A COMPUTATIONALLY EFFICIENT MODELLING OF 
LAMINAR SEPARATION BUBBLES 


FINAL REPORT 
for 

NASA Grant NAG-1-778 


Paolo Dini and Mark D. Maughmer 


Department of Aerospace Engineering 
The Pennsylvania State University 
233 Hammond Building 
University Park, PA 16802 


July 1990 



ABSTRACT 


In predicting the aerodynamic characteristics of airfoils operating at low 
Reynolds numbers, it is often important to account for the effects of laminar (tran- 
sitional) separation bubbles. Previous approaches to the modelling of this viscous 
phenomenon range from fast but sometimes unreliable empirical correlations for the 
length of the bubble and the associated increase in momentum thickness, to more 
accurate but significantly slower displacement-thickness iteration methods employ- 
ing inverse boundary-layer formulations in the separated regions. Since the penalty 
in computational time associated with the more general methods is unacceptable 
for airfoil design applications, use of an accurate yet computationally efficient model 
is highly desirable. To this end, a semi-empirical bubble model has been developed 
and incorporated into the Eppler and Somers airfoil design and analysis program. 
The generality and the efficiency have been achieved by successfully approximat- 
ing the local viscous/inviscid interaction, the transition location, and the turbulent 
reattachment process within the framework of an integral boundary-layer method. 
Comparisons of.the^predicted aerodynamic characteristics with experimental mea- 
surements for several airfoils show excellent and consistent agreement for Reynolds 
numbers from 2,000,000 down to 100,000. 
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Chapter 1 
INTRODUCTION 
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This thesis is concerned with a particular aspect of the overall task of estimating 
the aerodynamic performance of airfoils. Before a more precise statement of the 
problem can be made, it is necessary to describe its physical setting and to define 
some important terms. 

Background and Definitions 

The prediction of the aerodynamic performance of a two-dimensional airfoil 
can be obtained by letting the airfoil remain stationary while analyzing the flow of 
air over it. The speeds of interest here are usually sufficiently smaller than the local 
speed of sound of the air to justify an assumption of incompressible flow. 

As the Reynolds number, defined as 

r = (l.i) 

V 

falls below approximately R = 4 x 10 6 , the laminar boundary layer may run out 
of momentum before transition on the surface occurs. Since the momentum of the 
outer flow cannot readily reach the stagnant fluid near the surface, the imposed ad- 
verse pressure gradient can only be balanced by a negative velocity of the flow. That 

is, the boundary layer separates leaving a thin region of reversed flow underneath 

it. In such cases, transition occurs in the free shear layer downstream of laminar 
separation and is usually followed by reattachment as a turbulent boundary layer, 
such that a small amount of fluid remains trapped between the shear layer and the 
surface. This pocket of fluid is called a laminar separation bubble. 

Laminar separation bubbles on airfoils are observed over a large Reynolds num- 
ber range. Due to the ability of the turbulent shear layer to reattach to the surface 
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of a streamlined shape such as an airfoil, the critical Reynolds number is approx- 
imately two orders of magnitude smaller than that at which laminar separation is 
first observed. In the flow over a circular cylinder, a large jump in drag coefficient is 
observed as the Reynolds number is decreased approximately below R = 300,000. 
At this value, in fact, the bubble cannot reattach immediately downstream of tran- 
sition and massive laminar separation results. For airfoils, the critical Reynolds 
number depends on thickness and angle of attack but is usually between R — 40, 000 
and 200,000. Unlike the case of a cylinder, a bubble that has burst due to a de- 
crease in Reynolds number will usually reattach upstream of the trailing edge. In 
this state, it is called a long bubble and it causes a large decrease in lift and a 
large increase in drag. Bursting may also occur at high Reynolds numbers near the 
leading edge, downstream of the suction peak present on the upper surface at high 
angles of attack. The inability of the shear layer to reattach is ascribed in this case 
to the extreme values of the adverse pressure gradient. 

In the first forty years since laminar separation bubbles were first observed by 
Melville Jones in 1933, leading-edge bubbles at high Reynolds numbers received 
most of the attention since their bursting is generally believed to be responsible for 
abrupt stall. In more recent years, the development of small remotely controlled 
aircraft (Remotely Piloted Vehicles) has shifted the Reynolds number range of in- 
terest to values below one million. In this range, due to the delayed transition, 
bubbles may form near the mid-chord causing significant increases in profile drag, 
depending on their length and thickness. The shift from leading-edge to mid-chord 
bubbles has been a fortunate one. Since the latter are usually an order of magni- 
tude larger (although they appear to be very similar in internal flow structure), the 
greater ease of measuring pressure or flow variables has made it possible to make 
great progress toward full understanding of this phenomenon. This understanding, 
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however, is still not sufficient to establish a reliable bursting criterion. 

Whether bursting occurs due to a decrease in Reynolds number or to an in- 
crease in pressure gradient, the resulting flowfield is characterized by a global strong 
interaction between the viscous and the inviscid regions and section properties can- 
not be obtained by means of simple approximations. Away from these limiting 
conditions, however, the viscous/inviscid interaction induced by the bubble is lim- 
ited to its immediate vicinity, especially at low to moderate lift coefficients. Such 
bubbles have traditionally been called “short,” although it may be more accurate 
to refer to them as “weakly interacting.” 

Problem Statement 

The purpose of this study is to develop a model for weakly interacting laminar 
separation bubbles developing in the incompressible, two dimensional, viscous flow 
over airfoils that can accurately account for the increase in airfoil profile drag that 
accompanies them. 

From the point of view of boundary-layer theory, the distinction between weak 
and strong viscous/inviscid interaction is made on the basis of the magnitude of 
the modification to the inviscid pressure distribution by the viscous flow. A weakly 
interacting bubble, therefore, is such only insofar as its effect on the global invis- 
cid pressure distribution is concerned but is still characterized by a strong local 
interaction. In addition, when first confronted with the description of a region of 
recirculating flow between a separation and a reattachment point, it seems natural 
to assume that the flowfield is elliptic and that the boundary-layer approximations 
are not applicable. Given the very slow speeds of the recirculating flow and the 
thinness of the bubble, however, the boundary-layer approximations are usually as- 
sumed to be valid everywhere. Rather than within the recirculating flow, in fact, 
important signals are communicated upstream through the interaction with the 
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outer flow. A parabolic boundary-layer method is therefore adequate for analyzing 
laminar separation bubbles as long as it is coupled to the outer elliptic flow through 
at least a local interaction algorithm. This has been employed in the present model 
and, to maximize the computational efficiency, the boundary-layer development is 
calculated with an integral method. 

In order to determine the increase in drag caused by a bubble, the shear layer 
development must be calculated through the different regions of the bubble as shown 
in Fig. 1-1. The formation of a bubble is initiated at point <5, shown in the figure, 
by the laminar boundary layer separating from the airfoil surface. Using integral 
boundary-layer methods, this point can be determined with sufficient accuracy for 
airfoil design work. Once separated, the free shear layer development must be 
tracked and the transition from laminar to turbulent flow, which occurs near the 
point T, predicted. As shown in the figure, the separation bubble causes a plateau 
to form in the velocity distribution between the points corresponding to laminar sep- 
aration and the end of the transition region. From this point, the turbulent part of 
the bubble encompasses a pressure recovery region which leads to the reattachment 
of the turbulent shear layer at point K. As an additional pressure recovery always 
occurs downstream of a reattachment point [Green, 1966], the velocity distribution 
corresponding to the highly non-equilibrium, relaxing boundary layer downstream 
of reattachment usually “undershoots” the inviscid distribution. Eventually, the 
turbulent boundary layer reaches its fully developed state and the undershoot re- 
gion merges smoothly from below with the inviscid velocity distribution. Clearly it 
is possible, especially at low Reynolds numbers, that the turbulent boundary layer 
does not reach equilibrium before the trailing edge of the airfoil. 

Rather than describing the effect of the bubble on the pressure distribution, 
it is more informative to explain why such an effect is observed. This can be done 
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by invoking two conservation laws, the conservation of mechanical energy and the 
conservation of linear momentum. In an inviscid incompressible flow, the conserva- 
tion of mechanical energy is simply Bernoulli’s equation. While inside the boundary 
layer the total pressure is continually and unevenly being dissipated by the shear 
stresses, static pressure rise is still achieved by trading for it the available kinetic 
energy inside the boundary layer. Once the boundary layer runs out of kinetic 
energy shortly downstream of separation, since conservation of energy is a scalar 
law and energy is a positive-definite quantity, no more pressure can be recovered. 
Conservation of linear momentum, on the other hand, is a vector law such that, 
once the boundary layer runs out of momentum, it has no difficulty allowing it to 
become negative to balance the imposed inviscid adverse pressure gradient. As soon 
as this happens, however, the presence of the reverse flow, felt by the outer flow as 
a modification of the airfoil surface, effectively decreases this same inviscid pressure 
gradient. As the pressure gradient decreases, it, in turn, induces a lesser growth of 
reverse flow. In the limit, a bounded growth of reverse flow at constant pressure 
results. This new distribution of static pressure matches what the energy equation 
allows. 

While the momentum present in a laminar boundary layer is strictly dependent 
on what is “handed down” by the upstream development, the efficient cross-stream 
transfer of momentum by the Reynolds stresses that appear downstream of transi- 
tion in the bubble brings the momentum of the outer flow near the wall. This allows 
the reverse flow to be accelerated and near-inviscid pressure to be recovered by the 
reattachment point. This loss of outer flow momentum is observed, by definition, 
as a rapid boundary-layer growth in the turbulent part of the bubble. Within the 
context of an integral method, the momentum lost by the flow results in the large 
growth in momentum thickness measured in this part of the bubble and can result 
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in a significant increase in airfoil drag. 


Previous Models 


Previous attempts at modelling this flow phenomenon range from very crude 
empirical correlations for the transition length from laminar separation and subse- 
quent reattachment behavior to full simulations of the Reynolds-averaged turbulent 
Navier- Stokes equations. This variety of approaches is partly caused by the long 
time-span over which this problem has been studied. In the early years, several 
empirical models or empirical correlations for particular bubble characteristics were 
proposed. Von Doenhoff [1938] assumed that the dividing streamline is straight 
and tangent to the airfoil surface at the separation point. He determined the tran- 
sition location by assuming a constant transition Reynolds number, formed with 
the velocity at separation and the distance between separation and transition along 
the dividing streamline. The distance to reattachment is then determined using 
a constant spreading angle of the turbulent shear layer of 15° measured from the 
direction of the dividing streamline. Interestingly, the present model resembles this 
initial configuration more than many others that have followed in the next fifty 
years. Crabtree [1957] proposed a coefficient of pressure rise in the turbulent part 
of the bubble equal to 


Ptz — pr 

(1/2 )pV\ 


( 1 . 2 ) 


This coefficient has been used mainly to monitor bubble bursting, which would 
happen for values of a > 0.35. These early results, measurements, and hypotheses 
are discussed extensively in review papers by Ward [1963] and Tani [1964]. 

Gaster [1967] investigated bubble bursting with decreasing Reynolds number. 
He characterizes the bubble by the values of the momentum thickness Reynolds 
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number at separation, 


(R& 2 )s 


Usj^s 

V 


and a pressure gradient parameter, 


(82 )l A U 

v As 


(1.3) 


(1.4) 


where the velocity gradient in Gaster’s parameter refers to the mean inviscid gradi- 
ent between the separation and reattachment points. With varying conditions, for 
instance decreasing Reynolds number, a locus of points corresponding to the bubble 
evolution can be traced on a plot whose axes measure variations in these two param- 
eters. He found that bursting would occur always along the same line on this plot. 
Although bursting is not modelled in this study, these same two parameters play 
a fundamental role in the present model. Horton [1967] proposed a semi-empirical 
model where the governing integral boundary-layer equations are coarsely approxi- 
mated inside the bubble. This model will be discussed in more detail in Chapter 2. 
Van Ingen [1975] and Van Ingen et al. [1980, 1986] studied the bubble problem for 
many years. Their model is distinguished by a good approximation to the pressure 
distribution in the bubble region and of the transition process. The shear layer 
growth along the bubble, however, is not adequately approximated. These methods 
are usually not sufficiently accurate or general, their main shortcoming being an 
inability to account for the local ellipticity of the bubble flowfield. 

Since approximately 1970, advances in computer technology have prompted a 
number of more detailed viscous simulations by finite-difference methods. Briley 
and McDonald [1975] developed a local Navier-Stokes solution for the bubble region 
which is matched with the steady boundary-layer equations upstream and down- 
stream of the bubble and with the inviscid outer flow. Their predictions match 
Gault’s [1955] measurements to an acceptable degree of accuracy but comparisons 
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with airfoil drag data are not given. Davis and Carter [1984] developed an in- 
teracting method based on a perturbation of the outer inviscid flow by a source 
distribution representing the bubble and a finite-difference boundary-layer method 
for the bubble flowfield and boundary layer. By properly accounting for the local 
flow direction, their method is able to resolve rather interesting details about the 
flow inside the bubble. For instance, there appear to be three different vortices: 
one in the laminar part, one in the turbulent part, and one in between, next to the 
wall and rotating in the opposite sense. While the bubble characteristics predicted 
by this method compare well with the experimental data presented in their report 
[Gault, 1955], low Reynolds number cases and the effect of the bubble on airfoil 
drag are not considered. Cebeci and Schimke [1983], Cebeci [1989], and Kwon and 
Pletcher [1979] follow similar formulations, where a finite- difference boundary-layer 
method interacts with the outer flow. These methods are limited by the general- 
ity of the turbulence model they employ and, for the present time, are much too 
inefficient computationally to be used routinely. 

Over the same period of time, approaches that can be considered of an in- 
termediate degree of complexity have flourished: interactive methods that use an 
integral formulation for the boundary-layer development. Crimi and Reeves [1976], 
Drela and Giles [1987], Drela [1989], and Gleyzes et al. [1983] are of this type. The 
advantage of such methods is that they combine a relative computational efficiency 
with a potential for sufficient accuracy and generality. 

A Model for Airfoil Design 

In the Reynolds number range over which laminar separation bubbles form, 
40,000 < R < 4,000,000, they assume many different sizes and thicknesses, each 
tied to a particular airfoil geometry or to particular conditions on different airfoils. 
In fact, many different types of aircraft with different mission requirements and 
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design constraints fly in this range. It often happens that particular airplane design 
goals conflict with flow conditions that would minimize the detrimental effects of 
bubbles. For instance, the low-drag requirement of sailplanes necessitates extensive 
use of natural laminar flow technology. The favorable pressure distribution to 50 
or 60% of the chord generally employed on laminar flow airfoils necessarily leads 
to a steep pressure recovery over their aft-portion. The high performance of such 
airfoils hinges on forcing transition before the boundary layer encounters the ad- 
verse pressure gradient since the higher energy turbulent boundary layer is more 
able to recover the near-freestream trailing-edge pressure without separating. This 
is achieved by means of a region of moderately adverse pressure gradient, a tran- 
sition ramp,” at the end of the favorable pressure gradient region. As the Reynolds 
number decreases, the transition ramp becomes insufficient to destabilize the lam- 
inar boundary layer into transition before the beginning of the recovery such that 
a thick, high-drag bubble forms. The trade-off, therefore, is between the advan- 
tages of natural laminar flow at cruise speeds and the disadvantages of the bubble 
in thermalling flight. The higher c t required for thermalling would at first sight 
seem to help the designer in avoiding the bubble. In fact, higher ci s are obtained 
at higher angles of attack such that the adverse pressure gradient starts near the 
leading edge. The destabilizing effect of this type of pressure distribution is offset, 
however, by the lower Reynolds numbers characteristic of this flight regime. Since 
small changes in the pressure distribution can effect large changes in bubble struc- 
ture and, therefore, in drag, this “fine-tuning” engineering problem can be resolved 
only if the effects of the bubble can be accurately calculated under different types 
of pressure distributions at different Reynolds numbers. 

One way of accounting for laminar separation bubbles in airfoil design is the 
bubble analog used in the design and analysis program of Eppler and Somers [1980]. 
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The design method of this program uses an inverse conformal mapping method that 
allows great freedom and flexibility in specifying the characteristics of the airfoil 
pressure distribution to achieve the desired performance. The aerodynamic charac- 
teristics are calculated by an integral boundary-layer method driven by the inviscid 
velocity distribution. Close monitoring of the boundary-layer development is ac- 
tively used in the design process to aid in the modification of the inviscid velocity 
distribution to achieve the desired transition and separation behavior. In this pro- 
gram, the designer is warned about the presence of separation bubbles which might 
unacceptably increase the drag over that which is predicted assuming that transi- 
tion occurs at laminar separation. Although this approach has proven very useful in 
designing airfoils for low Reynolds number applications, it would be advantageous 
to have predictions of section properties which more fully account for the presence of 
laminar separation bubbles provided this can be done without significantly increas- 
ing computational requirements. In fact, while above R = 500,000 this criterion 
can be used effectively to design airfoils with short bubbles that do not increase 
the drag of the airfoil, as the Reynolds number decreases it becomes increasingly 
difficult to eliminate the detrimental drag increases [Mueller, 1984]. 

In order to design low Reynolds number airfoils more effectively, a new method 
of modelling the bubble has been developed in this thesis. It combines the speed of 
the empirical and semi-empirical approaches of the early years with the accuracy of 
the more recent interactive methods. This approach rests on the hypothesis that it 
should be possible to model a local phenomenon such as a laminar separation bubble 
through local rather than global information. This hypothesis has been confirmed. 
While the boundary-layer development upstream of laminar separation must in some 
cases be taken into account in order to predict accurately the transition location 
inside the bubble, all other characteristics of the bubble flowfield have been found 
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to depend on a few local scaling parameters that will be discussed in the following 
chapters. Furthermore, by accounting for the interaction through a local iteration, 
a uniformly accurate laminar separation bubble model has been developed. This 
model has been incorporated into the Eppler and Somers program. Although it is 
unable to account for strong interactions such as the large reduction in the suction 
peak sometimes caused by leading-edge bubbles, it is able to predict the increase 
in drag and the local alteration of the airfoil inviscid pressure distribution that are 
caused by bubbles occurring in the operational range which is of most interest. 

In Chapter 2, the reasoning leading to the independent parameters that control 
the bubble is retraced. In Chapter 3, the calculation of the laminar part of the 
bubble is described in detail. In Chapter 4, possible approaches to the modelling 
of transition are discussed together with the method employed here. In Chapter 5, 
the calculation of the turbulent part of the bubble is described in detail. Several 
empirical functions necessary to model this most complicated part of the airfoil 
flowfield are proposed. In Chapter 6, several airfoils are analyzed and the results 
are compared to experimental data. In Chapter 7, the range of validity of the present 
model is assessed, important results are summarized, and specific suggestions are 
given toward enlarging the empirical data base necessary to confirm the present 
formulation. 



13 


Chapter 2 

LOCAL PARAMETERS 


In this chapter, the approach followed in developing the present laminar sepa- 
ration bubble model is justified. It is then shown how the shortcomings of previous 
similar models can be remedied by letting the bubble flowfield depend on three local 
parameters. 


Modelling Philosophy 

The development of a model of a physical phenomenon entails a three-step pro- 
cess that is usually iterative rather than sequential: (1) identifying which dependent 
variables need to be modelled, (2) identifying which independent variables best char- 
acterize the conditions on which the phenomenon depends, and (3) determining the 
correct relationships between the two. In the case of incompressible fluid flow, (1) 
is comprised by the velocity and pressure field, (2) by spatial variables, time, and 
boundary conditions, and (3) by a differential relationship, the Navier-Stokes equa- 
tions, which embodies pointwise mass and momentum conservation. This model can 
be integrated numerically to obtain the dependence of the velocity and pressure on 
varying geometry and flow conditions for the laminar separation bubble problem as 
well as for thousands of other flowfields. 

Motivated by computational efficiency requirements, methods of varying de- 
grees of approximation have been developed. As soon as approximations are intro- 
duced, the modelling challenge changes in nature since assumptions have to be made 
about what does and what does not need to be approximated and these assumptions 
must be supported either by analytical arguments or by experimental evidence. For 
instance, viscous/inviscid interaction methods rely on the boundary-layer assump- 
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tions and on the extensive analytical and experimental evidence that supports them. 
There are two important consequences of this type of approximation: (1) the num- 
ber and types of flowfields that can be analyzed with interactive boundary-layer 
methods is severely restricted from what could be done by the original equations, 
(2) by integrating the inviscid or outer problem, the independent conditions are 
shifted from geometry and spatial variables to flow variables. The final solution is 
arrived at by allowing the outer and boundary-layer flows to interchange their roles 
as dependent and independent processes at each iteration, with the relationship 
between them given by the numerical integration of either Laplace s equation or 
the boundary-layer equations. Although the solution information is limited to the 
airfoil surface, approximating high Reynolds number viscous flows in this way adds 
valuable integral insight about the solution over and above the knowledge that it is 
governed, pointwise , by the Navier-Stokes equations. 

If the computational requirements are even more stringent, as in the present 
case, then a more drastic approximation of the physical phenomenon is necessary. 
The resulting model will be even more limited in applicability and it will have to 
rely on relationships that are further removed from the original differential pointwise 
balance and closer to the integrated solutions. In essence, the limit to this approx- 
imation process is simply an explicit solution to the original mathematical model 
of the phenomenon, necessarily obtained at the expense of its generality. That is, 
if such a solution is not possible at a particular level of approximation, the model 
is simplified further. The advantages of speed and insight brought by an explicit 
solution, however, are off-set by a model that may have become too simplistic and 
restrictive in applicability. In this thesis, the approximation is brought one step be- 
yond the viscous /inviscid interaction approach to a semi-empirical method. In this 
case, this approach has proven to be a successful compromise between speed, gener- 
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ality, and understanding. In order to minimize the reliance on empirical guesses, 
the governing integral boundary-layer equations are enforced in the bubble and are 
complemented with a set of simple relationships that capture the essential features 
of the bubble allowing the prediction of its characteristics and of its effects in gen- 
eral. This formulation was arrived at by starting with the simpler models and 
adding complexity only after establishing that it was absolutely necessary. 


Early Results 

In the course of the research reported here, efforts to develop a method able to 
predict the effects of a laminar separation bubble which interacts weakly with the 
inviscid flow began with the incorporation of the classical empirical model of Horton 
[1967], modified according to the suggestions of Roberts [1980] and Schmidt and 
Mueller [1989], into the Eppler and Somers program. Because they are formulated 
in terms of integral boundary-layer properties, bubble models such as these are well 
suited to the integral boundary-layer analysis method employed by Eppler. This 
method employs two coupled governing differential equations, the momentum and 
energy integral equations, 


d 6 2 
ds 


£/ _ (it .O') — — 

2 ^ Hl 2 +T> U ds 


d&z _ „ Jj_dU 
ds Cd 6 U ds 


( 2 . 1 ) 

( 2 . 2 ) 


together with appropriate closure relations for c/, Cd, and H \2 [Eppler, 1963], 
given in the Appendix. Contrary to simpler, one-equation methods such as that 
of Thwaites [1949], in two-equation methods the shape factor (# 32 * in this case) is 
obtained directly from the governing equations and is therefore independent of the 
local pressure gradient parameter. This allows such methods to analyze accurately 
the non-similar boundary-layer developments characteristic of aerodynamic flows 
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provided that the assumed family of velocity profiles approximates the actual flow 
reasonably well. 

The empirical models noted above do not take advantage of the accuracy af- 
forded by a two-equation method to calculate the development of the shear lajer 
along the bubble. Instead, they obtain the growth of S 2 along the bubble by means 
of rather coarse approximations of these equations. Thus, assuming a constant- 
pressure plateau between separation and transition and negligible skin friction 
brought by the near-stagnant fluid in this region leads to 

(<$2)r = (^2)5 (2-3) 

from the momentum integral equation. Based on low Reynolds number measure- 
ments [Schmidt and Mueller 1986; O’Meara, 1986; Brendel, 1988], Schmidt and 
Mueller [1989] suggest using, instead, 



from the similarity solution for the laminar free shear layer. The value of momen- 
tum thickness growth predicted by this equation, which was proposed earlier by 
Russell [1978], increases with decreasing Reynolds number. In order to evaluate 
this equation, the length of the laminar part of the bubble, given by correlations to 
be discussed in Chapter 4, must be known. 

From the transition point, the growth of S 2 in the turbulent part of the bub- 
ble is approximated by simplifying the momentum and energy integral equations. 
Combining these two equations leads to Truckenbrodt’s shape parameter equation, 

6 2 ^ = (H 12 - 1 )H 37 S -L^.+C D -fx 3 2 (2.5) 

As discussed by Horton [1967], it appears from experimental data that 



dE 3 2 
ds 
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Using this result along with that of vanishing skin friction at a point of reat- 
tachment, Horton was able to reduce Eq. ( 2 . 5 ) to 


'62 dU' 


°d 1 

Uds 

71 

[h 32 (H 12 -1)J 


= Arc 


( 2 . 7 ) 


The assumption of a universal reattachment velocity profile and a constant Cd in the 
turbulent part of the bubble leads to a constant value for A-jz. Assuming constant 
Cr> and H-i'i and a linear pressure recovery from transition to reattachment, one 
can integrate the energy integral equation to obtain the momentum thickness at 
reattachment, 


(62)11 = ($2 )r 





( 2 . 8 ) 


Then, eliminating (^2)7^ between Eqs. ( 2 . 7 ) and ( 2 . 8 ) one obtains 


~ (life) 1 - (fe) 4 


( 2 . 9 ) 


In order to implement this result numerically, Utz is decreased in small increments 
from the value of Ur- At each step, £ 2 is calculated and it is checked whether or 
not the segment joining T io'R. intersects the inviscid velocity distribution. When 
this happens, (62)21 is calculated from Eq. ( 2 . 8 ) and the turbulent boundary-layer 
method is started at s-ji using this value and (££32)11 = 1-51 for initial conditions. 

This type of formulation is inadequate for several reasons. It does not properly 
account for the effects of the local viscous/inviscid interaction on the pressure dis- 
tribution in the laminar part and therefore cannot accurately calculate the growth 
of 62 in this region. It relies on local empirical transition criteria which are un- 
able to sense the influence of the upstream boundary-layer development. Also, in 
actuality, neither H32 nor Cr> are constant in the turbulent part of the bubble. 
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The turbulent recovery distribution can deviate appreciably from a straight line. 
Finally, the pressure at the reattachment point can be significantly below or above 
the corresponding inviscid value. 

In an attempt at better approximating the pressure distribution induced by 
the bubble, van Ingen’s functions were implemented [Ingen and Boermans, 1986]. 
In the laminar part use is made of the empirical function 

— = .978 + .022exp(— 4.545£ - 2.5£ 2 ) (2.10) 

Us 


where 




S — Ss 

(R 6i )s(S 2 ) s 


( 2 . 11 ) 


which allows for a slight pressure recovery between separation and transition. For 
the turbulent part, a Stratford pressure distribution is employed and again it is 
assumed that reattachment occurs at the intersection with the inviscid distribution. 
Fig. 2-1 shows a comparison between the Horton and Stratford recoveries. The 
dashed line is the locus of possible reattachment points as given by Eq. (2.9). 

Using the empirical separation bubble model noted, the sensitivity of the 
boundary-layer development and drag prediction to various parts of the bubble 
was explored. In order to achieve accurate drag polar predictions, it was found 
necessary to capture the vanishing of the bubble with increasing adverse pressure 
gradient while at a constant chord Reynolds number. It is found that the transition 
length for such empirical models responds to variations in chord Reynolds num- 
bers but not in pressure gradient. The drag prediction, in turn, is very sensitive 
to small variations in the governing parameters, for instance, the pressure level at 
the beginning of the turbulent pressure recovery. Thus, although generally capa- 
ble of predicting features of the bubble to within thirty percent, empirical bubble 
models based only on conditions at separation cannot provide acceptable drag pre- 
dictions. Fig. 2-2 shows the aerodynamic characteristics of the Eppler E387 airfoil 
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at R = 300, 000 obtained with this early version of the bubble model compared to 
the original Eppler and Somers prediction and the measurements of McGhee et al. 
[1988] taken in the Low Turbulence Pressure Tunnel at the NASA-Langley Research 
Center. 

Before beginning a discussion of possible alternatives to the above approxima- 
tions, two issues indirectly affecting the bubble model should be addressed. Firstly, 
due to the presence of the boundary layer, the experimental pressure distribution 
does not equal the inviscid one at the same angle of attack but falls inside it, leading 
to a smaller lift coefficient. If a weakly interacting bubble is present, it will modify 
the viscous pressure distribution only locally. Using the inviscid pressure distribu- 
tion to drive the bubble model, therefore, necessarily leads to a discrepancy between 
the predicted and measured results if the same angle of attack is prescribed. This 
discrepancy can be eliminated by employing a viscous/inviscid interaction algo- 
rithm. While appropriate for analyzing near-stall conditions, methods of this type 
are not really necessary at lower lift coefficients where the design effort is most often 
concentrated. In fact, since aerodynamic characteristics are usually compared at 
the same ct rather than at the same oq the difference in angle of attack between the 
experimental and inviscid lift coefficients poses no obstacles to comparing drag pre- 
dictions obtained with the inviscid pressure distribution to the experimental drag 
polar. It may be expected that the validity of the present model will decrease for 
mid-chord bubbles on highly aft-loaded airfoils and for leading-edge bubbles near 
bursting. 

Secondly, the turbulent boundary-layer analysis method developed by Eppler 
and employed in the program is based on empirical equilibrium relationships be- 
tween the integral variables [Eppler, 1963]. While quite appropriate for analyzing 
high Reynolds number flows without bubbles, such a method cannot correctly ac- 
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count for the relaxing, nonequilibrium turbulent boundary layer downstream of 
reattachment, especially at lower Reynolds numbers. The present model, therefore, 
makes use of the nonequilibrium turbulent boundary-layer method developed by 
Drela [1986]. To provide a comparison between these two boundary-layer meth- 
ods, in Figs. 2-3 and 2-4 the Eppler E3S7 and the NASA NLF(1)-1015 airfoils 
are analyzed assuming transition at the laminar separation point using the origi- 
nal and Drela’s turbulent boundary-layer methods. The turbulent separation point 
predicted by Drela is downstream of that predicted by Eppler for the low Reynolds 
number case but is equal to it at the higher Reynolds number. The drag coefficient, 
however, is consistently higher by 10 to 15 counts. 

The Local Independent Parameters Controlling the Bubb le 
It is stated in Chapter 1 that the present model has confirmed the hypothesis 
that it should be possible to model the weakly interactive bubble solely by means of 
local parameters, with the exception of the transition process. Although transition 
may be the most important effect, it is not of much help without an accurate 
estimation of the bubble flowfield. Thus, the failure of the early empirical models 
can be traced not only to their poor modelling of transition but also to their inability 
to capture the principal physical processes in the two parts of the bubble. In spite of 
this, their reliance on local conditions has not been abandoned in the present model. 
Rather, it has been extended to achieve a more detailed prediction of this flowfield. 
Specifically, the local parameters must be able to characterize three aspects of the 
bubble flowfield that have been found to control all other bubble characteristics: 
the boundary-layer momentum thickness at separation, the mean inviscid pressure 
gradient along the bubble, and the thickness of the bubble at transition. In addition, 
the correct lowest-order response of the bubble to these inputs must be identified. 
The failure or success of previous approaches to predict the bubble flowfield and its 
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Fig. 2-3 Comparison of the Eppler and Drela turbulent boundary-layer methods at a low 
Reynolds number. Transition at laminar separation. 
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effects on airfoil performance can be shown to be directly related to their failure or 
success to account properly for these independent and local flow conditions. 

The first investigators of the bubble problem recognized that the transition 
location plays a fundamentally important role in determining the size of the bubble. 
They therefore recognized the need to estimate the degree of instability present in 
the boundary layer as it separates from the surface. A very stable boundary layer 
at separation, in fact, is likely to be followed by a fairly long bubble, whereas an 
unstable one would be associated with a shorter bubble. The Reynolds number 
has been universally used as an indicator of how far from turbulence a particular 
flowfield is. In addition, as the chord Reynolds number is decreased, the bubble 
becomes longer until it bursts at the critical value. It seemed plausible, therefore, 
that a good correlating parameter for the stability of the separating boundary-layer 
flow could be arrived at by forming a Reynolds number with the value of inviscid 
velocity and momentum thickness at laminar separation as characteristic velocity 
and length, 

(Rs 2 )s = — ^ (2.12) 

This parameter is certainly useful but, unfortunately, only provides a coarse indica- 
tion of the stability of the flow at laminar separation. Another parameter that has 
been used extensively in empirical correlations is the value of momentum thickness 
at laminar separation nondimensionalized with respect to the airfoil chord. Rather 
than a measure of the stability of the separating shear layer, this parameter only 
contains information on how much momentum has already been lost by the bound- 
ary layer upon reaching the laminar separation point. These criteria lack sensitivity 
to the effect of the pressure distribution on transition, namely the effect of the up- 
stream boundary-layer development. As will be discussed in Chapter 4, this effect 
is well modelled by the e n method of linear stability theory. 



26 


Another effect that needs to be accounted for is the local strong viscous/inviscid 
interaction. Since this interaction is a result of the presence of separated flow and 
since the amount of separated flow is affected by the severity of the pressure gradient 
and by (^ 2 ) 5 , it is not surprising that this effect has been found to scale well with 
a parameter that incorporates these flow conditions, Gaster’s pressure gradient 
parameter, given by Eq. (1.4). 

It should be mentioned at this point that a recent study [Pauley et al., 1989] 
of the laminar separation bubble shows Gaster’s parameter to be important in yet 
another respect. In this study, the unsteady laminar Navier-Stokes equations are 
discretized to calculate the flow through a duct. Although the local Reynolds num- 
ber of the flow based on the distance between the entrance of the tunnel and the 
laminar separation point is on the order of R = 50, 000-300, 000, transition to turbu- 
lence cannot be accurately computed (computational unsteadiness can occur) due 
to the grid resolution and time-step used to keep the computational requirements 
within reasonable limits. A suction port on the upper wall of the duct provides 
an adverse pressure gradient for the laminar boundary layer developing along the 
lower wall. 

The value of the parameter P mai = *24, a modification to Gaster’s parameter 
with the velocity gradient representing the maximum rather than the average in- 
viscid gradient between separation and reattachment, is found to correspond to the 
boundary between steady and unsteady reattachment of the laminar shear layer. 
This boundary correlates well to Gaster’s bursting line such that it separates the 
long (steady, P ma x < .24) from the short (unsteady, P max > -24) bubbles measured 
by Gaster. Because of this correlation between the steady and long bubbles and the 
unsteady and short bubbles, it is proposed by Pauley et al. that the reattachment 
process may be governed by the large-scale laminar pressure field rather than by 
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turbulent transfer of momentum. While this cannot be proved or disproved with 
certainty, the onset of unsteadiness can be justified in the case of laminar flow. In 
fact, as Pmax increases, the inviscid gradient and/or the momentum loss at sepa- 
ration increases. Thus, a greater momentum transfer is necessary to recover the 
inviscid pressure and/or to accelerate the sizable amount of reverse flow. As the 
laminar shear layer is not capable of such momentum transfer, reattachment prob- 
ably occurs due to Coanda’s effect, or the formation of a low-pressure region below 
the shear layer. This effect as it relates to laminar separation bubbles is discussed 
by Russell [1978]. When the adverse inviscid gradient exceeds the favorable suction 
of the flow below the shear layer, the bubble starts growing without bounds; until, 
that is, the large vortex of recirculating flow at constant pressure becomes unstable 
and is shed causing the bubble to collapse in size. This small bubble then starts 
growing again and the cycle repeats. 

Before the above description was arrived at, it was thought that transitional 
bubbles, too, were unsteady in the large scale and short-time mean. It was also 
thought that the unsteadiness near reattachment would feed back upstream and 
influence the transition location. It seems now that if there is some unsteadiness 
it should not necessarily be periodic in nature and that the turbulent transfer of 
momentum is a much more significant factor in determining reattachment than 
Coanda’s effect if turbulence is present. Although the final word on transition can- 
not be given with the method used in the present model, it appears at this point to 
depend mainly on upstream rather than downstream conditions. It is still not clear, 
therefore, why the unsteady laminar simulation reported by Pauley et al. correlates 
so well with G aster’s measurements of transitional bubbles. Resolution of this issue 
may have to wait for a complete understanding of bubble bursting. 

Given a means of approximating the viscous/inviscid interaction in the laminar 
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part of the bubble and the transition location, the height of the bubble at transition 
can be estimated. This geometrical characteristic of the bubble is the keystone that 
bridges the laminar and the turbulent parts, providing the correct characteristic 
length for the latter. As the spreading angle of the turbulent shear layer is a 
weak function of Reynolds number and pressure gradient, in fact, the length of the 
turbulent part of the bubble follows directly and the model is thereby closed. 
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Chapter 3 

THE LAMINAR PART OF THE BUBBLE 


The Eppler and Somers program uses a very reliable criterion to detect laminar 
separation. It is based on the value of the energy to momentum thickness shape 
factor, 

(tf 32 )s = 1.515095 (3.1) 

This value is approached from above. Upon detection of laminar separation, the 
development of the separated laminar shear layer is calculated using the momentum 
and energy integral equations, together with closure relations to be discussed below. 
Instead of implementing this boundary-layer method in the inverse mode as it is 
usually done, the development of a general family of pressure distributions in the 
laminar part of the bubble allows its calculation in the direct mode. 

Pressure Distribution 

The function used to approximate the pressure distribution in the laminar part 
of the bubble is a generalization of that developed by van Ingen and Boermans [1986] 
and given by Eq. (2.10). This distribution allows a slight pressure recovery after 
laminar separation, quickly approaching a limiting value. Using detailed pressure 
distributions in the bubble region available from wind-tunnel tests of the NASA 
NLF(1)-1015 airfoil in the NASA-Langley Low- Turbulence Pressure Tunnel, the 
accuracy of Eq. (2.10) was checked for several different conditions. It was found 
that, as the pressure gradient along the bubble decreases, the pressure distribution 
tends to fall below van Ingen’s curve while, as the pressure gradient steepens, it 
becomes flatter, closer to Horton’s approximation and above van Ingen’s curve. It 
is therefore postulated that Eq. (2.10) can be improved by relaxing the amount of 
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pressure recovery between separation and transition, 

— = (1 - DU) + DU exp(— 4.545£ - 2.5?) ( 3 - 2 ) 

Us 

The steeper the pressure gradient along the bubble, the smaller the value of DU . 
This behavior is consistent with an inviscid velocity distribution calculated over an 
ever- thickening displacement surface in a steepening adverse gradient. 

While the agreement with the experimental distributions was much improved 
by use of Eq. (3.2), an inconsistency became apparent when attempting to predict 
the pressure distribution over leading-edge bubbles. Given the very large gradients 
along these bubbles, the predicted pressure distributions in the laminar part were 
quite flat, in contrast to the measurements which show only a small perturbation 
of the inviscid distribution with a significant pressure recovery between separation 
and transition. This apparent contradiction with the trend observed for mid-chord 
bubbles can be resolved once it is realized that the amount of pressure recovered 
is inversely proportional to the magnitude of the perturbation of the displacement 
surface, that is, to the amount of fluid entrained by the bubble. Near the leading 
edge, the boundary layer is so thin that the short bubbles (1 — 6%c) usually observed 
there can only hold a very small amount of fluid and can therefore only modify the 
inviscid distribution slightly. 

Although there is a strong correlation between the thickness of the boundary 
layer at separation and the amount of fluid caught in the separated region, the 
value of (<^ 2)5 more precisely reflects the input to the momentum balance that 
determines such amount: the greater the momentum already lost by the boundary 
layer, that is, the greater the amount of separated flow necessary to counteract the 
imposed inviscid gradient. The variable DU should therefore depend both on the 
average pressure gradient along the bubble as well as on (£ 2 ) 5 - Both these effects 
are included in G aster’s pressure gradient parameter, Eq. (1.4) which, for this 
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reason, is thought to be a better independent parameter than simply the average 
dimensionless inviscid velocity gradient along the bubble, 


AjU/Uoc) 

A(s/c) 


(3-3) 


In dimensionless variables P becomes, 


P = R 


(62). 


l 2 


A{U/Uqq) 

A(s/c) 


(3.4) 


From experimental pressure distributions, it is found that DU is well repre- 
sented as a function of the G aster pressure gradient parameter, P. This functional 
relationship, shown in Fig. 3-1, was developed by extracting corresponding values 
of DU and P directly from the experimental pressure distributions of the NLF(l)- 
1015 [NASA LaRC LTPT, June 1987] and the Eppler E387 airfoils [McGhee et ah, 
1988]. The solid line is a quadratic least-squares fit that has been included in the 
model, 


DU = 


0.0610 + 0.3048P + 0.5072P 2 -P < .3 
0.0152 -P > -3 


(3.5) 


The value of DU = 0.022 used by van Ingen and Boermans falls in the middle of 


the variation in DU shown in Fig. 3-1. 

As pointed out by van Ingen [1989], the factor 4.545 in Eq. (3.2) was derived to 
ensure continuity in the velocity gradient when Thwaites’s laminar boundary-layer 
method is used to determine the laminar separation point. In fact, if the derivative 


of Eq. (2.10) is evaluated at laminar separation and the variables are rearranged, 
one obtains 

= — (4.545)(.022) = -.10 (3.6) 

5 

which is Thwaites’s laminar separation criterion. Eppler’s separation criterion (Eq. 


'(S 2 ) 2 dU 

v ds 


(3.1)), however, corresponds to slightly different values of Thwaites’s parameter 
for differing upstream developments. Furthermore, in the present formulation a 
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variable DU is used in place of 0.022. Therefore, the constant necessary to ensure 
a continuous velocity gradient at separation when Eq. (3.2) is used in conjunction 
with Eppler’s laminar boundary-layer method is 


C = 


When C is substituted for 4.545 in 


DU 


'(* 2 ) 2 dir 

v ds 


Eq. (3.2), the result is 


(3.7) 


— = 1 — DU { 1 — exp 
Us l 


' 1 U's 

DU U s 


( s - 



(3.8) 


where the prime denotes the derivative with respect to 5 and the second term in 
the exponent of the original expression has been neglected. 


Closure 

In order to integrate the momentum and energy integral equations in the direct 
mode as driven by the pressure distribution given above, closure correlations for 
H 12 , Cf , and Cd must be provided. The most natural choice is to use the reversed 
Falkner-Skan, or Stewartson [1954], profiles since the attached Falkner-Skan, or 
Hartree, profiles [Schlichting, 1979] are used to develop the correlations upstream 
of separation. Recent measurements by Fitzgerald and Mueller [1990], however, 
seem to indicate that the Stewartson profiles may not be the best choice. This 
matter was therefore examined in some detail. 

Fitzgerald and Mueller [1990] have obtained good agreement between their 
measurements and the two-parameter profile family originally developed by Green 
[1966] for a turbulent shear layer forming a free stagnation point downstream of 
a base. As shown in Fig. 3-2, the two parameters are linked to the geometrical 
characteristics of the profiles. ( h/b ) is the ratio of the distance of the shear layer 
from the centerline of the wake to the width of the shear layer and G is the am- 
plitude of Coles’s wake function. Since there is slip along the centerline of such a 



Fig. 3-2 
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recirculating base flow, these profiles cannot be used to develop a correlation for c/. 
By applying the definitions for the integral thicknesses of the boundary layer and 
for the dissipation coefficient, the following relationships are obtained: 


(l-§G)-2f(l-2G) 

(2 - f G + f G 2 ) + (1 - 3G + 2 G 2 ) 

(1-|G)-2|(1-2G) 


(3.9) 

(3.10) 


R> 1 C D =^f-[l-la + 2- b (l-2G)] (3.11) 


In order to compare these relationships to those obtained from the Stewartson 
profiles, it is necessary to know how the two parameters vary inside the bubble. The 
values used by Fitzgerald and Mueller to fit the profiles measured inside one bubble 
can serve as a starting point. The three boundary-layer variables are evaluated 
at values of the parameters corresponding to the same downstream station inside 
the bubble and then plotted against one another. The same calculations are then 
repeated for values of G and (h/b) similar to those used by Fitzgerald and Mueller 
in order to determine the sensitivity of the correlations to these parameters. The 
result is shown in Figs. 3-3 and 3-4 where these new two-parameter correlations 
are compared to those developed by Drela [1986] from the Stewartson profiles. The 
solid lines utilize the fitted variations of G and (h/b). As both H\2 and H 32 increase 
monotonically between separation and transition, moving to greater values of the 
abscissa on these plots corresponds to moving downstream inside the bubble. Thus, 
both are very similar to the Stewartson correlations near separation but can be 
quite different further downstream. 

While # 12 (^ 32 ) seems quite sensitive to changes in the parameters, 
Cd{H 32 , Rs 7 ) is not, thereby making the determination of its exact dependence 
on G and (h/b) less crucial. It appears from the measurements that the back- flow, 


GREEN PROFILES 
FITTED 



Comparison of the shape-factor correlation for the Falkner-Skan and Green profiles 
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proportional to G , may be constant within each bubble although different for differ- 
ent bubbles. As shown in Fig. 3-3, the values of shape factors actually measured, 
although different in absolute value, follow the same slope, thus confirming a con- 
stant value of back-flow velocity. These considerations justify eliminating ( h/b ) 
between Eqs. (3.9) and (3.10) and expressing the closure relationships in terms of 
Hz 2 i calculated from the governing equations and G, whose behavior within each 
bubble appears easier to correlate to local flow conditions, 


H 12 = 


R(,.Cd 


3(1 — G) — Hz2 
(l-G)(l-2 G) 


7T 


! G 3 




(4-5G)(l-G)-(2-3G)ff 3 ; 
4(1 - G) - 2# 3 2 


(3.12) 

(3.13) 


Although very encouraging, these results do not appear sufficiently well devel- 
oped to be implemented in the model in their present form. In fact, the dependence 
of G on local flow conditions is unknown and no measure of c/ can be obtained 
from these profiles. By contrast, although the details of the velocity profiles are 
not well represented by the Stewartson profiles, the integrated parameters derived 
from them are not far from the corresponding values obtained with the fitted Green 
profiles. It seems better, therefore, to keep using the Stewartson profiles, for now, 
until these issues have been resolved. Accordingly, the closure correlations devel- 
oped by Drela [1986] have been implemented in the model. Since the governing 
integral equations yield the value of H ^2 directly, Drela’s shape factor correlation is 
inverted, 

'H 32 -1.194068^ 2 


TT #32 - 1.194068 , 

Hu = SSI + 




.04 


V 


-64.4 


(3.14) 


The skin-friction coefficient is found from 

-0.067 + 0.01977 # 12 < 7.4 

-0.067 + 0.022 [l - , #12 > 7.4 




£/ _ 

2 


(3.15) 
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while the dissipation coefficient is given by 

R 6a — = 0.207- 0.003(i?i2 - 4) 2 (3.16) 

H 32 

Removal of the Separation Singularity 

Knowledge of the pressure distribution downstream of separation has led to a 
simple technique for removing the Goldstein singularity at laminar separation and 
to allow for some effect of the bubble on the pressure distribution upstream of lam- 
inar separation. When the development of the laminar shear layer is calculated in 
the direct mode from laminar separation using the function described above, the 
boundary-layer variables remain at their separation values for a few percent chord 
before starting to grow normally. This is believed a consequence of the singularity, 
in the following sense. The growth of the boundary-layer variables upstream of sep- 
aration does not reflect just the local pressure gradient but is increasingly affected 
by the singularity as separation is approached. This effect is equivalent to a much 
steeper pressure gradient than is actually present and causes the distribution of H 32 
to exhibit a very steep slope immediately upstream of separation. The skin-friction 
coefficient behaves similarly and thereby leads to a prediction of the separation 
point upstream of the experimentally observed location. Therefore, when the sep- 
aration values are incremented downstream of separation using a pressure gradient 
continuous with its value upstream, it is felt by the boundary layer as, in fact, a 
much gentler gradient which cannot maintain the previous rate of growth. In the 
present model, the point where the separated laminar shear layer starts growing 
again is taken as the actual laminar separation point. The laminar boundary layer 
is therefore calculated again from a few percent chord upstream of the “inviscid” 
laminar separation point to this point by prescribing between them an assumed 
(cubic) development of H\ 2 {$) and solving the laminar boundary-layer equations 
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in the inverse mode. Smooth growth of the boundary layer through the separation 
point results. 

The local inverse solution causes the velocity distribution to start deviating 
from the inviscid some distance upstream of the laminar separation point. As this 
distribution matches the experimental measurements, this behavior of the separa- 
tion point and of the corresponding velocity distribution is believed to be correct. 
In fact, as the boundary-layer assumptions break down as separation is approached, 
a modification to the inviscid pressure distribution upstream of separation is to be 
expected. Referring to Fig. 3-5, as the new separation point is usually at a higher 
pressure than the original, the old value of DU is usually too great and a new one 
is calculated from 

DU new = DUoid - (3.17) 

( Us)old 

Using this value and the new value for the velocity gradient at separation, Eq. (3.8) 
can be used to generate a new velocity distribution whose tangent is continuous 
with the current distribution at the separation point. 


Separation Angle 

Using the new value of momentum thickness at the “viscous” separation point, 
the tangent of the angle the separating streamline makes with the surface is given 
by an empirical relationship proposed by Wortmann [1974], 


tan 7 = — 




(3.18) 


A similar relationship proposed by van Ingen et al. [1980] was also examined, 


tan 7 = 


( Rs 2 )s 


(3.19) 


B is given the experimental mean value of 17.5. The model did not perform at all 
well with the latter relationship. The reason is that the factor 64P assumes a wide 




Close-up of the laminar separation region demonstrating the removal of the 
Goldstein singularity. 
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range of values, from as low as 2 up to 30 for the range of bubbles examined, and 
their average does not give the model enough flexibility. Eq. (3.18), on the other 
hand, gives the correct scaling for this variable. In addition, very similar values 
for 13 as given by Eq. (3.18) have been reported by Pauley et al. [1988] in their 
Navier-Stokes simulation. Eq. (3.18), therefore, has been included into the model. 


Chapter 4 
TRANSITION 
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The prediction of the transition location inside the bubble has received a great 
deal of attention since the very first models. In fact, both the drag increment due to 
the bubble as well as the bursting behavior depend strongly on the location where 
the calculations switch in a more or less gradual way from laminar to turbulent 
flow. Although it has been claimed that the length of the transition region inside 
the bubble must be modelled accurately [Walker, 1989], a point-transition has been 
found to work very well in the present model. In this chapter, a few empirical 
transition criteria are discussed together with the method employed in the model. 


Empirical Criteria 

Ever since bubbles were first observed, many researchers have looked for an 
empirical correlation between the distance from separation to transition and local 
bubble characteristics such as conditions at separation. Since it was observed that 
the length of the bubble is inversely proportional to the Reynolds number, the first 
transition criterion, proposed by Von Doenhoff [1938], assumed a constant Reynolds 
number based on the velocity at separation and the distance between separation 
and transition, 

R t = = 50,000 (4.1) 

V 

As the Reynolds number increases, Us usually increases, too, such that this criterion 
forces the laminar length to decrease in size. Horton [1967], thirty years later, used 
a value of 40, 000. This expression can be rearranged as 


h 


40,000 


(62)3 


c 


c 


(4.2) 
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O’Meara, and Mueller [1986] took an experimental mean from many sets of data as 

(4.3) 


* = 155 <^ 
c c 


which is a less general relationship. As pointed out by Schmidt and Mueller [1989], 
these relationships cannot capture the vanishing of the bubble at the experimentally 
reported value [Crabtree, 1957] of momentum thickness Reynolds number at sep- 
aration of about 750, over which transition precedes laminar separation. Schmidt 
and Mueller describe how the correlation developed by Vincent de Paul [1972] is 
able to capture this effect by allowing a variable Similarly, based on the same 
data sets used by O’Meara and Mueller, Schmidt and Mueller propose the following 
criterion. 


ft f P13- 3.7820 (Re, )s]^, 27 < (R s ,) s < 72 

C \ [ 267 - 0 . 3709 (^) 4 ]^, 72 < (R t , )s < 720 ' 1 

such that the bubble length vanishes at values of (R$ 2 )s > 720. Fig. 4-1 shows 

this criterion in graphical form together with Eq. (4.3). Since Eq. (4.3) does not 

incorporate any dependence on (Rs 2 )s , it is shown in the figure as a plane whose 

projection onto the (8 2 )s / c-Z\ / c plane is simply a straight line through the origin. 

While it is true that the Reynolds number has a strong influence on the tran- 
sition length, and therefore on the overall length, of separation bubbles, this effect 
is reflected in many of these correlations in a rather coarse way. The correlations 
for which the transition length is proportional to the value of momentum thickness 


at separation better predict the length of leading-edge bubbles due to a spurious 
coincidence. For airfoil flows, in fact, small values of (# 2)5 correspond to leading- 
edge bubbles, which form on large suction peaks near the stagnation point and are 
usually quite short. Larger values of (S 2 )s are usually associated with mid-chord 
bubbles that occur far downstream of the stagnation point and are usually much 
longer. Fig. 4-2 illustrates this effect, where the isolated points are leading-edge 
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Fig. 4-1 Comparison of two empirically derived transition correlations 
i for the separated laminar shear layer. 
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+ O’Meara and Mueller 
O Schmidt 


Fig. 4-2 Comparison of the laminar lengths predicted by the two 
correlations for the Eppler E387 airfoil at R = 300,000. 
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bubbles. The value of (<*>2)s> however, does not in itself contain any information 
about the stability of the separating shear layer. In fact, no empirical criterion 
developed so far can capture the effects of a variable pressure distribution. Specif- 
ically, no local criterion can capture the vanishing of the upper-surface mid-chord 
bubble with increasing angle of attack. 

Eppler’s Transition Criterion 

Given the success of Eppler’s transition criterion in predicting the transition 
location in attached boundary layers, it was thought that it should be possible to ex- 
tend such a criterion to separated boundary layers. Whereas the criteria described 
above seek a correlation between transition in the free shear layer and conditions 
at separation, which is at a different location on the airfoil, Eppler’s criterion is 
based on the local characteristics of the boundary layer. Since, unlike in the early 
bubble models, the boundary layer development here is calculated also downstream 
of separation, it seemed that monitoring the development at each downstream in- 
crement would naturally lead to a more accurate transition prediction. In order 
to better explain how Eppler’s criterion is implemented, Eppler’s boundary-layer 
development plot should be described. 

Integration of the momentum and energy integral equations gives the values of 
82 and 83 at each downstream station. The shape factor H32 = 83/82 is therefore 
also known. Since the inviscid velocity along the airfoil is taken as the boundary- 
layer edge velocity that drives the boundary-layer development, U and 82 at each 
downstream station can be grouped to obtain the development of R & 2 . Eppler con- 
nects subsequent (H32 , i ?£ 2 )-pairs on a plot whose axes measure the variation in 
these two variables, thereby describing the boundary-layer development from the 
stagnation point to the trailing edge in a very concise way. The stagnation point 
occurs at a value of H32 = 1.62 and Rs 2 approximately equal to 10 . Values of 
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Rs 2 < 10 are not plotted. Fig. 4-3 shows a typical boundary-layer development. 
Following the upper surface development given by the solid line, laminar separa- 
tion is encountered when H Z 2 = 1.515095. This criterion is shown on the plot as a 
vertical dotted line. From the laminar separation point, H 32 starts growing again 
until transition is met. While P32 grows monotonically inside the laminar part of 
the bubble with downstream distance from laminar separation, R $ 2 stays approxi- 
mately constant. A very similar criterion to Eq. (4.4) can therefore be constructed 
by correlating the value of #32 at transition to (Rs 2 )s- The particular transition 
criterion for separated shear layers shown in this figure as a family of cubics intro- 
duces an additional dependence on P . Thus, while the bubble is seen to disappear 
at values of (Rs 2 )s > 875, for smaller values (#32 )t increases with decreasing P. 
This trend in (P32 )r does not necessarily imply a similar trend in laminar length 
since the rate of growth of #32 depends on the pressure distribution in the lami- 
nar part of the bubble. Although unsuccessful, this and similar criteria served to 
illustrate how poorly (#32 )r correlates to local bubble conditions incorporated in 
parameters such as P and ( Rs 2 )s * In fact, while any one particular bubble could be 
matched quite easily with a criterion as shown in the figure, any such criterion was 
consistently found of very limited generality. Downstream of transition, the tur- 
bulent shear layer growth causes H Z 2 to decrease again to a value at reattachment 
which is weakly dependent on the local momentum thickness Reynolds number and 
always close to 1.51. The attached turbulent boundary-layer development causes 
Hz 2 to increase again toward the flat-plate fully developed value of 1.77. As the 
trailing edge is approached, the adverse pressure gradient drives the boundary layer 
toward turbulent separation which occurs at # 32 = 1.46 for Eppler’s method and 
#32 = 1 .51 for Drela’s method, which will be discussed in the next chapter. 

In the course of these investigations, it was found that the curve corresponding 
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Fig. 4-3 Modified Eppler boundary-layer development plot showing the shear la 
ment inside the bubble. 
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to attached transition will cause natural transition to be predicted in many cases 
where a bubble is actually present. This happens because the value of Reynolds 
number at which transition precedes laminar separation can be quite high. In fact, 
it can be higher than the value of 720 proposed by Schmidt and Mueller. Eppler’s 
original transition criterion [Eppler, 1963], not shown, intersected the laminar sep- 
aration line at R& 2 = 463. Based on additional measurements, Eppler [1989] has 
modified the criterion to 

[In R St } r > -21.74 + 1SAH 32 + 125 {H i2 - 1.573) 2 (4.5) 

where now the intersection occurs at R$ 2 = 704. For such a value, however, tran- 
sition is still predicted too soon when the boundary layer is near separation. In 
order to be able to compare bubbles forming on airfoils at chord Reynolds numbers 
of one million or greater, approximately, the 125 in Eq. (4.5) is replaced with 190. 
This leads to a value at the intersection of R$ 2 — 875, which seems to work well. 
This apparent coarseness in Eppler’s criterion is due to its having been calibrated 
mainly from airfoil drag coefficient data rather than from a modelling of the actual 
transition process. The reason why transition is predicted instead of a bubble at 
high chord Reynolds numbers is simply that at such values the bubble does not 
cause any significant drag increase; in fact, it usually causes the same increase in 
momentum thickness that is given by the attached turbulent boundary layer over 
its length. 

In an attempt to reconcile Eppler’s transition criterion with another very suc- 
cessful transition prediction method, the e n method, they are plotted together in 
Fig. 4-4. In the e n method, a value of n = 9 has been found to correspond 
to observed transition locations in many different flows. Rather than using Ep- 
pler’s boundary-layer development plot, this comparison is done here on a plot of 
Rs 2 vs. H 12 in order to include several measurements of these variables at observed 
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|Tjg 4-4 Comparison of Eppler’s transition criterion with the amplification surface for self- 
similar profiles. Dashed line for H 12 < 4 is Eppler s transition criterion for attached 
flow; dashed line for H \ 2 > 4 is one of the attempted criteria for separated flow. 
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transition locations. This comparison was attempted also because, since the e n 
method can predict transition successfully in a separated shear layer, it was hoped 
that a suitable criterion for use in the boundary-layer development plot could be 
inferred from it. In this figure, contours of constant n are given for the Falkner- 
Skan self-similar developments [Schlichting, 1979]. For such developments, in fact, 
a unique surface exists that allows the determination of the value of n from the local 
values of the momentum thickness Reynolds number and shape factor alone. The 
equation defining this surface was developed by Drela [1986] and will be discussed 
below. Shown in this figure is Eppler’s transition criterion for attached boundary 
layers as well as one of the attempted criteria for separated flow. It is interesting 
that Eppler’s curve falls quite close to the n = 9 contour, for zero-pressure gradient 
flow (i?i 2 = 2.59). As can be seen from the wide scatter in the experimental data, 
however, it seems dubious that a single curve could capture the correct transition 
location inside the bubble with any generality. In fact, the growth of n in a non- 
similar boundary layer development will not follow the surface whose contours are 
shown in the figure, regardless of whether the flow is attached or separated. This 
implies that transition is not dependent solely on the local boundary-layer charac- 
teristics but depends also on the manner in which the boundary layer arrives at such 
values. The inclusion of path-dependency, or the effect of upstream boundaxy-layer 
development on transition, has captured the generality that the transition criteria 
discussed so far lack In order to fully understand why, the e n method will now be 
discussed in detail. 


The e n Method 

The e n semi-empirical transition prediction method, developed thirty years 
ago [van Ingen, 1956; Smith and Gamberoni, 1956] and since successfully applied 
to a variety of aerodynamic flows, relies on a theory that can explain the moderate 
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success of some of the correlations described above, that can correctly distinguish 
between leading-edge and mid-chord bubbles, and that can correctly model the ef- 
fects of variations in pressure distribution upstream of the bubble. Linear stability 
theory, in fact, directly models the growth of instabilities in a boundary layer while 
indirectly, through the boundary-layer development, accounting for the effects of 
Reynolds number. n(s) is defined as the logarithm of the ratio of distuibance am- 
plitude at station s to its amplitude at neutral stability, s 0 . Transition is assumed to 
take place when n reaches a value previously correlated to experimentally observed 
transition locations. For similar flow environments this value has been reported 
to lie around 9 by many researchers, although it appears to depend on Reynolds 
number [Evangelista and Vemuru, 1989; Horstmann et ah, 1990]. 

It is generally accepted that linear stability theory correctly models the transi- 
tion process for approximately 70% of the distance between neutral stability (n = 0) 
and fully turbulent flow. The actual “transition region,” however, is usually defined 
as the region between the first appearance of turbulent spots and fully turbulent 
flow [Arnal, 1984], or the last 30% of this distance. In order to approximate the 
transition process in the nonlinear amplification region, an intermittency function 
is usually employed at a value close to n = 8. Equally often, transition is taken to 
be completed when n — 9 — 14 and the turbulent calculations are started abruptly 
at the corresponding streamwise station. In both cases, significant empirical input 
is necessary to render the method usable. Another concern associated with this 
method is that the amplitude of the disturbance at neutral stability is unknown. 
The inability of the method to model directly the influence of the turbulence in- 
tensity of the oncoming air and surface conditions on the onset of transition forces 
the use of further empirical corrections [Mack, 1977]. In spite of these weaknesses, 
the e n method remains the engineering transition prediction method most faithful 
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to the actual physical process. 

The e n method is based on the numerical solution of the Orr-Sommerfeld 
equation which is derived from the Navier- Stokes equations as follows. The two- 
dimensional Navier-Stokes equations are perturbed and linearized. An assumption 
of locally parallel mean flow is made. Subtracting the mean flow, a system of three 
partial differential equations for the perturbation field is obtained. Since the coef- 
ficients are functions only of y, the method of normal modes can be applied in the 
particular form 


v'(x,y,t) = Re 1 " 


(4.0) 


which implicitly assumes a sinusoidal disturbance. Here v l is the disturbance in the 
y-direction. Similar expressions are assumed for u f and p } . a* is the wavenumber 
and to* the radian frequency, which can be nondimensionalized with respect to a 
reference length and velocity, 


oc = a*L re f 

oj*L ref (4-7) 

CO = — 

Uref 

Substitution of these expressions into the partial differential equations for the dis- 
turbance field results in a set of three ordinary differential equations for the dis- 
turbances. Eliminating u' and p' in favor of v ' , the Orr-Sommerfeld equation is 
obtained, 


v"" + [~iR(aU - u) - 2a 2 ]v" + [iR(aU - u)a 2 + iRU"a + a 4 ]u = 0 (4.8) 


or and u are, in general, both complex. Although the derivation and the analysis 
are performed in the complex plane, it is understood that physical quantities are 
obtained by taking the real part of any of the complex variables employed. 

Two limiting cases of interest are 

f a = a r , to — oj r T iw,- temporal instability , . 

a = a r + iai , oj — u> r spatial instability ' ’ ' 
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To calculate the amplification of disturbances as the boundary layer develops, a 
spatial instability analysis is necessary. Thus, the Reynolds number and the fre- 
quency are specified and the wavenumber is found by solving numerically the Orr- 
Sommerfeld equation subject to the boundary conditions 
, dv' 

v = -j— = 0 at y = 0 and as y — * oo (4-10) 

The homogeneous boundary conditions imply the existence of an infinity of solutions 
or eigenvalues a. The value of greatest physical interest is the eigenvalue for which 
— a,- is largest. This can be seen from the expression for the fluctuation, 

v' = v(y)e i[(Q *+' Q *) r -^ t ] 


[v(y)e 


~ Q i *1 ,,«(“* *-"r0 


(4.11) 


where v(y ) is the distribution of the disturbance amplitude (eigenfunction) and 
the bracketed term is the amplitude of the disturbance at a distance x from some 
reference point. The ratio between this term evaluated at two different ^-locations 
represents the amplification of a disturbance between them. Assuming a boundary 
layer developing in zero pressure gradient, 


v(y)e 


-Of; X 2 


M = 

Ax u(y)e - °'i II 


_ (x 2 -*l) 


(4.12) 


The amplification factor, n, is defined as the logarithm of this ratio, 


J = -a*(x 2 -xx) (4.13) 

The amplification factor is thus equal to the area under the amplification rate curve. 
In this case, this curve is a constant because the parallel flow approximation and 
the absence of a pressure gradient prevent any change in the local Reynolds number 
such as, for instance, R^ 2 . 

The continuous wavelength modulation encountered by a fixed-frequency dis- 
turbance travelling downstream in a developing boundary layer is approximated 
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locally by a series of constant-wavelength plateaus, obtained by solving the Orr- 
Sommerfeld equation for a sequence of parallel mean flows of different thicknesses 
and at different pressure gradients, each characterizable by a different value of i?$ 2 
and i?i 2 - As the streamwise extent of the plateaus tends to zero, Eq. (4.13) can be 
generalized to 

n = / —a* dx (4. 

Jxi 

As a measure of how far from transition a boundary layer is, it is preferred to 
calculate the total amplification that has occurred: the amplification, that is, from 
the lower branch of the neutral curve, 



n 



—a* ds 


(4.15) 


where the independent variable is again s, the distance along the airfoil from the 
front stagnation point. 

In general, the disturbance environment in a wind tunnel or in free flight does 
not consist of a single-frequency pressure pulse but, rather, of broad-band or white 
noise. The boundary layer simply amplifies those frequencies that, at a specific 
local Reynolds number, correspond to what in simpler dynamical systems is called 
the natural frequency. Since the local boundary-layer Reynolds number changes 
with downstream distance, different frequencies are amplified as the boundary layer 
develops. This necessitates the tracking of several different frequencies at the same 
time, with the first that reaches n = 9 indicating transition. 

As the frequencies amplified in the separated shear layer may be different from 
those upstream of separation, it may not be necessary to monitor the stability of 
the laminar boundary layer upstream of separation. The fact that some frequencies 
may come close to n — 9 upstream of separation does not necessarily imply that 
transition will occur sooner once the boundary layer separates. If this is true, 
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then it should be possible to devise a transition criterion based solely on local 
information, for instance on conditions at laminar separation. The consistent failure 
of all previous such criteria would still not be sufficient proof of the untenability of 
this hypothesis. 

Prelaws Approximate e n Method 

While the questions raised above need further study, for the present time the 
approximation to the e n method developed by Drela [1986] has been implemented 
in the Eppler and Somers program. In order to demonstrate that this method intro- 
duces an error in the calculation of n over and above its declared approximations, 
it is now discussed in some detail. 

Rather than performing a linear stability analysis of the boundary-layer velocity 
distribution as can be obtained, for instance, from a finite-difference method at each 
downstream station, following Gleyzes et al. [1983] Drela computes a data base of 
the stability characteristics of the Falkner-Skan profiles that can be “tapped” during 
a boundary-layer calculation using the local shape factor as the coupling parameter. 
More precisely, the nondimensional growth rate, —a,*, corresponding to a particular 
value of the local shape factor of a Falkner-Skan profile and of the local Reynolds 
number is divided by the local boundary-layer characteristic thickness, i.e. 62 , 
to obtain the physical growth rate to be used in the integral (4.15). Given that 
the correct characteristic thickness is obtained independently, from the momentum 
integral equation, the manner in which the nondimensional data base is generated 
is of no consequence. The most convenient way is to calculate the growth rates for 
self-similar developments at constant H 12 values and increasing Rf> 2 . Starting from 
the Orr-Sommerfeld spatial instability analysis of the Falkner-Skan profiles at many 
different values of H\ 2 , a set of neutral curves is generated, one for each value of 
shape factor. An example of these neutral curves is shown in Fig. 4-5 for two values 
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of H i 2 - For each value of shape factor, the dimensionless amplification rate, —a,-, 
is evaluated along rays of constant reduced frequency, 

2nfv 


F = 


U 2 


(4.16) 


to form curves — ati(Hi 2 , Rs 2 ,F). From these curves the amplification factor for the 
development of each Falkner-Skan profile is found from 
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where m is constant and defined as 
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Denoting the momentum-thickness integral (not a function of s ) by I, 
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Thus, Equation (4.17) becomes, 

n(Hi2,Rs 2 ,F) — — f -(a*S 2 ) dRt 2 
•' R '> 20 

= 7777^ I**’ dR S2 (4.22) 

17 ( 7712 )] Jr 62q 

In this way, the original dimensionless eigenvalues of the dimensionless Orr- 
Sommerfeld equation can be used to find n, which is defined in terms of a di- 
mensional wave number and distance. This can be done for a self-similar profile 
since dRs 2 /ds assumes the particular form shown. 

The curves obtained with this integral for different frequencies and at a constant 
H \2 are shown in Fig. 4-6. Drela takes the envelope as a straight line as done by 
Gleyzes et al. [1983]. This leads to the following expression for the amplification 
surface for self-similar developments, 
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where the superscript “e” denotes a value obtained from the envelope of amplified 
frequencies and 
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While Gleyzes et al. then evaluate the amplification integral as 
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Drela changes the variable of integration back to s, 


n(s) = f 

j s 0 


dn 

dR$. 


(Hu) 
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(4.27) 


At this point, Drela finds an expression for dR$ 2 /ds in a rather roundabout way. 
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He shows that dRs^/ds equals this last expression with a 1 in place of the sec- 
ond term inside the square bracket. If this is taken as a condition, the resulting 
differential equation can be integrated as follows, 
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(4.29) 

Thus, the functional form characteristic of Falkner-Skan profiles is recovered, where 
the constant is in general a function of the shape factor. Drela now introduces the 
two “empirical relations” (given in the Appendix), 
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What he means is “analytical correlations from the Falkner-Skan profiles.” Since 
he is using these profiles here and since his form for dRs 2 /ds indeed implies the 
assumption of a Falkner-Skan ^-growth, it is not clear why he did not simply 
substitute the Falkner-Skan expression for <5 2 (s) directly into his expression for 
dRs 2 /ds , Eq. (4.28), 
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which is identical to Eq. (4.21). Thus, Drela only had to curve-fit /(fiT^). 

Having established that 

+ = (4.33) 


Drela’s integral for evaluating n(s), 
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can be written, using Eqs. (4.17), (4.21), (4.26), and (4.33), 
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where 8 2 (s) comes from the non- similar boundary-layer development as calculated 
by the momentum and energy integral equations. Thus, 
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Which is identical to the expression used by Stock and Degenhart except for the use 
of an envelope to find or,*. Using as calculated by the governing equations is 

precisely what enables the method to account for upstream history on the growth of 
n. That is, the dimensionless growth rate —a* obtained at each downstream station 
from the value of H 12 and R& 2 is divided by the local boundary-layer momentum 
thickness which is in general different from the value in a self-similar development 
at the same value of i7 a2 and R$ 7 , 

The Envelope Error 

In collaboration with Selig [1990], the author found that, even if Drela had used 
the actual envelope of the amplification curves instead of approximating it with a 
straight line, his n(s)-development would not correspond to the true envelope of 
the amplification curves at constant frequency along the airfoil surface. The error 
arises whenever the boundary-layer development is non-similar. To see this, it is 
necessary to rely on a numerical example, since the functions in question are not 
known analytically but are obtained numerically. Thus, it is helpful to envision a 
fictitious boundary-layer development made up of two constant- H\ 2 lengths with a 
discontinuous jump in between. Fig. 4-5 shows the neutral curves corresponding to 
the two values of shape factor. It is desired to compare the growth of n obtained 
by following the development through the jump in H 12 at constant frequency to 
that obtained using Drela’s envelope. Fig. 4-6 shows the amplification curves for 
the three reduced frequencies shown on the neutral curves plot as calculated by Eq. 
(4.22) together with the envelopes given by Eqs. (4.23)-(4.25) for the two values 
of H\ 2 . Fig. 4-7 shows the n-growth along the boundary layer with the switch in 
shape factor occurring at R$ 2 ~ 500. 

The three frequencies selected represent limiting cases that serve best to eluci- 
date the argument. Referring to Fig, 4-7, as R$ 2 increases n(Fi) grows according 





Amplification curves for two values of the shape factor. 
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Fig. 4-7 Growth of n in a non-similar boundary-layer development compared to Drela’s 
erroneous envelope method. 



66 


to Fig. 4-6 up to the maximum and, just as it is ready to start decaying, the jump 
in H \2 forces further amplification until the upper branch of the neutral curve cor- 
responding to H \2 = 2.67 is crossed. This additional growth will not necessarily 
be steeper than the envelope. n(F 2 ) does not start being amplified until the switch 
occurs, at which point it grows quite steeply in accordance with the greater area 
under the amplification rate surface — a'i(2.67, i?$ 2 ,u;). This curve does not nec- 
essarily exceed the envelope. Starting with the F$ -curve, at all lower frequencies 
the growth of n will follow the H \ 2 = 2.67 line, which is parallel to but greater 
than Drela’s envelope. In this example, it is possible to recover the steep similarity 
growth given by Eq. (4.23) since the shape factor is held constant downstream of 
the switch. In a non-similar development, however, the variation of H\ 2 is con- 
tinuous. If a monotonically increasing shape factor is approximated by a series of 
infinitesimally small steps, the resulting growth on n will never be able to “catch 
up” with the value obtained from a self-similar profile at the same local shape factor 
and Reynolds number. The correct envelope obtained by following each frequency, 
therefore, will lie above Drela’s approximation without ever reaching the growth 
given by Eq. (4.23). The converse is true for an accelerating boundary-layer. 

Based on the above argument, Drela’s envelope method may be expected to 
overpredict the transition location for non-similar, decelerating flows and to under- 
predict it for non-similar, accelerating flows. This method has been incorporated 
in the model nonetheless, although the method of Stock and Degenhart should be 
extended to separated profiles and used in its stead in future research. 
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Chapter 5 

THE TURBULENT PART OF THE BUBBLE 


The calculation of the turbulent part of the bubble relies on the assumption that 
reattachment will occur. An independent bursting criterion has not been devised, 
nor have existing ones been tested. The reason is that bursting occurs either at 
very low Reynolds numbers or when the mean inviscid pressure gradient becomes 
too steep, downstream of a suction peak. Regardless of the fact that the same 
mechanism may not be responsible for both types of bursting, both conditions 
represent extremes that lie outside the capabilities of the simple approach taken with 
the present model. This is not so much because of a failure of the bubble model itself 
but, rather, because at such extremes the onset of strong global viscous/inviscid 
interaction modifies too greatly the inviscid pressure distribution which drives the 
model. 

In this chapter, a parameter that characterizes the turbulent part of the bubble 
is introduced and the modelling of the reattachment process within the context of 
an integral method is discussed in detail. 

Scaling Parameter 

Having obtained a good approximation of the strong local viscous/inviscid 
interaction induced by the laminar part of the bubble and a fairly accurate transition 
location, now a steep pressure recovery must be predicted in order for the turbulent 
shear layer to reach again the inviscid distribution as it flows past reattachment, 
downstream of the strong interaction region. In the laminar part, the strength of the 
interaction depends on the amount of near-stagnant fluid downstream of separation 
and can be gauged by the deviation of the local viscous pressure distribution from 
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the inviscid. As this amount is, in turn, proportional to the local mean inviscid 
pressure gradient and to the momentum already lost by the boundary layer — to 
(£2)5? that is — it is not too surprising that the pressure recovered in this part of 
the bubble correlates well with P. In the turbulent part, on the other hand, the 
strong interaction is the result of a different mechanism, which acts to reverse what 
happened in the laminar part. As a consequence, the solution may be expected to 
depend on a different scaling parameter. 

Assuming the boundary-layer equations to be valid in this region allows the 
use of conditions at transition as initial conditions for the turbulent calculations. A 
similar approach to that taken in the laminar part was attempted, at first. Thus, 
several types of velocity distributions, similar to Horton’s straight line or to Strat- 
ford’s recovery, were prescribed and the boundary-layer equations were solved in the 
direct mode. In addition to the difficulty in approximating observed reattachment 
velocity distributions with any degree of generality, the solution was found highly 
sensitive to the smallest variations in the input pressure distribution, suggesting 
that this part of the bubble could not be calculated by solving the boundary-layer 
equations in the direct mode. 

Unlike in the laminar part, it is more convenient here to approximate the 
distribution of #32 than that of edge velocity. In fact, given that the value as well 
as the slope of the #32 -distribution is always known at the reattachment point, 
a general function has been developed in this study which allows the solution of 
the turbulent part of the bubble in the inverse mode. The distribution of #32 is 
specified as 


#32 (y) — (#32)72 + [(-#32)7" ~ (#32)72.] 


j" 1 + sin 



( 5 . 1 ) 


where the subscript i = 1,2 denotes the amplitudes of the sin(l/x) function up- 
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stream and downstream of reattachment, respectively, and 


y = Q - yoj cr + y 0 

_ 7T 

37r — sin _1 (l/>li — 1) 

f (s - St )/£ 2 S < STZ 

\ [(5 — sj )/£ 2 — ljSFd- 1 s > s-jz 


(5.2) 

(5.3) 

(5.4) 


where 


SF = \fMjM 


(5.5) 


ensures continuity in the curvature of # 32 ( 3 ) at the reattachment point. This 
function is shown in normalized form in Fig. 5-1. 

The inverse boundary-layer formulation employed here, where the distribution 
of shape factor is specified [Eppler, 1989], is especially convenient and powerful since 
it allows complete control of the boundary-layer behavior in an otherwise extremely 
sensitive region while at the same time relying on an intrinsically general function. 
Indeed, in view of the discussion given above concerning the reversal of solution 
hierarchy in regions of reversed flow from the standard weakly-interacting boundary- 
layer formulation, it is perhaps not mere coincidence that an inverse method should 
prove so much more effective in this part of the bubble. Such effectiveness, however, 
comes at a price. In fact, while specification of the pressure recovery distribution, 
if the correct one were indeed known in general, would automatically drive the 
boundary layer to reattach at the correct location, the turbulent length of the 
bubble, £ 2 -, in Eq. (5.4) is not known a priori and must be found by independent 
means. 

For some time during the model development, the inviscid velocity distribution 
was used to guide the location of the reattachment point. The correct value of £ 2 , 
that is, would be the one that leads the pressure at reattachment near the inviscid 
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value at the same station. Previous models, such as Horton’s, have also assumed 
that reattachment occurs at the intersection of the bubble recovery with the inviscid 
pressure distribution. Using the experimental pressure distributions for the NLF(l)- 
1015 airfoil, however, integration of the boundary-layer equations in the direct mode 
along the bubble revealed that the reattachment point may lie significantly above 
as well as below this intersection. Although the validity of this result is limited by 
the accuracy of the closure correlations employed, as well as by the assumption of 
the validity of the boundary-layer equations themselves, the observed trends seemed 
too consistent to be mere coincidence. Specifically, it was observed that for long, 
high-drag, mid-chord bubbles, corresponding to the middle of the airfoil drag polar, 
the reattachment point is always below the inviscid distribution (above, in the value 
of pressure) whereas for short mid-chord bubbles about to disappear corresponding 
to the top of the low-drag bucket in the airfoil polar, as well as for leading-edge 
bubbles, the reattachment point is always well above the inviscid. Allowing this 
point to move and trying to reproduce the observed trends by means of an empirical 
function based on local bubble conditions was helpful in obtaining reasonable drag 
predictions at least for a single Reynolds number. The inability to justify such a 
function physically and its limited generality, however, prompted further study. 

It was noticed that interaction methods have no trouble pinpointing the unique 
solution. Whereas the accuracy of such computed solutions still depends on the ac- 
curacy of the correlations, their mere ability to reach a converged solution indicates 
that all the necessary physical constraints are somehow accounted for. Also, their 
success at predicting the correct trends with respect to the location and pressure 
level of the reattachment point is explained by their ability to “sense” the presence 
of the wall through variations in the strength of the transpiration velocity or through 
the shape of the displacement thickness itself needed for convergence. An equiv- 
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alent geometrical constraint had to be introduced in the present model. Whereas 
prescribing the pressure level at reattachment is consistent with a direct formula- 
tion of the boundary-layer equations, in a region of such strong viscous/inviscid 
interaction as the bubble the inverse formulation of the boundary-layer equations 
should be complemented with the treatment of the correct physical process as the 
independent one to reflect the reversal in the solution hierarchy. The local flow- 
field is driven by a turbulent momentum-transfer mechanism whereby the outer 
momentum brought toward the wall accelerates the near-stagnant reverse flow un- 
til reattachment is achieved. Reattachment, therefore, becomes dependent on the 
efficiency of this mechanism. In geometrical terms, the reattachment location be- 
comes dependent on the spreading angle of the turbulent shear layer and on the 
initial distance of the shear layer from the wall; that is, on the height of the bubble 
at transition. Since the spreading rate of the shear layer is nearly insensitive to 
variations in Reynolds number and, in this particular flowfield, in inviscid pressure 
gradient, the turbulent length of the bubble becomes almost entirely dependent on 
the thickness of the bubble at transition nondimensionalized with respect to the 
airfoil chord. Making the characteristics of the turbulent part of the bubble depend 
on this parameter has enabled the model to reproduce all the available experimental 
data with excellent and consistent accuracy in the range 100, 000 < R < 2 , 000, 000. 


Governing Equations 

Following Eppler [1989], the distribution of #32 described above is input into 
the momentum and energy integral equations expressed in the inverse mode, 
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where denotes the derivative with respect to s of Eq. (5.1). 

Drela’s turbulent boundary-layer method makes use of an additional equation, 
the turbulence lag equation. Following Green et al. [1973], Drela simplifies the 
stress-transport equation originally proposed by Bradshaw and Ferris [1968] to a 
rate equation for the maximum shear stress coefficient, 


6 (1C T 
C T ds 



(5.8) 


where C T is defined as 


C r = 


( — u'v') r 


and 6 , the boundary-layer thickness, is given by 


(5.9) 


8 — 62 ^3.15 + — — ^ + $1 (5.10) 

This system of equations needs to be supplemented with closure relations for Co , 
# 12 , c/, and CV e ,. 


Closure 

Much of the following discussion is an expansion on that given by Drela [1986]. 
The maximum shear stress coefficient is used in the expression for the local dissi- 
pation coefficient. Its deviation from the equilibrium value, C Tcq , as calculated by 
means of Eq. (5.8) accounts for the slow response of the intensity of the turbulence 
being convected from upstream to varying local conditions. The equilibrium value 
refers to the equilibrium turbulent boundary layers of Clauser [1954], for which the 
local pressure gradient acting on the displacement thickness is balanced by the local 
wall shear stress. The ratio of the two forces acting on an incremental “slice” of the 
boundary layer is a constant, 

<$i dp 2 5\ dU 

t w ds Cf U ds 


(5.11) 
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For equilibrium flows, the modified shape parameter Gt is also constant, 

_ #12 ~ 1 1 1 82 

' Hn s/^n y/^J/2 A ( • 1 

A is the displacement thickness scaled in an analogous way to the defect profile, 

f°° U - 
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(5.14) 


To lowest order, the defect law (u — f7)/u* vs. y/8 will collapse the outer layer 
of any boundary-layer development for which equals a constant. Thus, there 
exists a unique relationship between Gt and (3 t for equilibrium boundary layers, or 
“equilibrium locus” [Kline et al., 1968], 

Gt = 6.7VlTa75^ (5.15) 

Using Eqs. (5.11), (5.12), and (5.15), the velocity gradient can be expressed as 

(5.16) 
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Noting that for these flows H 12 and H $2 are nearly constant, Eq. (2.5) can be 
simplified to 

62 dU 
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Eliminating the velocity gradient between Eqs. (5.16) and (5.17), yields an expres- 
sion for the dissipation coefficient valid only for equilibrium flows, 
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(5.18) 


In the laminar boundary-layer calculation, the closure relationships are derived 
from the similarity profiles and, strictly, are valid only for such flows. This means 
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that if for a self-similar development the Falkner-Skan pressure gradient parameter 
is calculated directly from its definition [Schlichting, 1979], 

v ds 

this value will be constant and will correspond to a unique value of the shape factor, 
also constant. In a non-similar development, on the other hand, the relationship 
between j3 and H 32 is not unique anymore. It is found in this case that calculating 
the shape factor from the governing equations while disregarding entirely the local 
value of (3 allows accurate values for Cd , H 12 , and c/ to be obtained even though 
the similarity correlations are utilized. In the calculation of the turbulent boundary 
layer the same approach is used by Eppler [1963] who utilizes empirical correlations 
that are valid mostly for equilibrium flows. In this case, however, obtaining the 
shape factor from the governing equations is not sufficient anymore if the pressure 
gradient is varying too rapidly as, for instance, downstream of a bubble. This is 
because of the large inertia of the Reynolds stresses that respond slowly to variations 
in local pressure gradient. As the stress level is the most important physical quantity 
in a turbulent boundary layer, it makes sense for there to be a need for its careful 
modelling in flows that depart too greatly from equilibrium assumptions. Just 
as in the laminar case decoupling the shape factor from local conditions brought 
a great gain in accuracy, so here the maximum shear stress coefficient should be 
decoupled from a statement of the type of Eq. (5.21) (in the next page) and obtained 
independently. Therefore, just as in the laminar case an additional equation is 
necessary, the kinetic energy integral equation, a third governing equation must be 
introduced, the rate equation for C T . 

The integral variable which is affected the most by the stress level in the bound- 
ary layer is the dissipation coefficient. In the laminar case, obtaining the shape 
factor independently was sufficient to obtain accurate values of Cd through the 


(5.19) 
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similarity correlation, Eq. (3.16) inside the bubble. Although this provides no 
guarantee of success for the turbulent case, an analogous argument is followed and 
Eq. (5.18) is assumed to be valid also for non-equilibrium flows. In order to be able 
to use it, however, it must be expressed in terms of C Teq - To this end, the dissipation 
coefficient is assumed to equal the sum of a wall and a wake contribution, 

C D = c f U 9lip + 2C r . f (1 - U a ii p ) (5.20) 


By equating this expression to Eq. (5.18), 
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So that, finally, the dissipation coefficient is given by 


Cd = CfUslip + 2C r (l — Ualip ) 


(5.21) 

(5.22) 


(5.23) 


where U a u p is given by Eq. (5.22) and C T is calculated from the rate equation (5.8). 

The shape factor correlation is derived from the analytical profiles of Swaf- 
ford [1983]. In the present method, this correlation needs to be expressed in the 
form 77 , 2 ( 7732 ). Since Drela’s correlation is expressed as 7732(77,2) and, unlike in 
the laminar case, cannot be inverted, a close approximation has been developed. 
Defining first 


Hz 2 0 — 1*505 + 

62 

(5.24) 

400 

77 , 2 0 — 3 + 

Tt5 2 

(5.25) 

c, = 0.081(7?*, -300) 0 1 

(5.26) 

c 2 = 0.0158(7?*, -300) 0 08 

(5.27) 

3000 

C3 -l-06 + ( +666)1-5 

(5.28) 



77 


H 12 is obtained from 
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This expression, as the original, is valid for > 400. It is shown in Fig. 5-2. 
The skin-friction law of Swafford is employed for attached flow, 
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(5.30) 


This expression does not give the correct value of skin-friction in the turbulent part 
of the bubble. The modification employed is discussed in the next section together 
with the remaining empirical functions. 


Supplementary Functions 

In order to achieve good agreement with measured bubble pressure distribu- 
tions and corresponding airfoil drag coefficients, existing turbulent correlations had 
to be modified and several empirical functions have been introduced. These func- 
tions are specified by assigning values to certain parameters left free. While such a 
formulation granted the model great flexibility, the unknown dependence of these 
parameters on local flow conditions significantly limited its generality. Identification 
of the correct scaling parameter for the turbulent part of the bubble, however, has 
made the determination of the correct functional dependence of these parameters 
both easier and less crucial. The remaining details for the calculation of the tur- 
bulent part of the bubble will now be discussed and, wherever possible, supporting 
arguments will be given. It should be realized, however, that final proof of the 
validity of the modifications introduced will await more detailed measurements of 
this region of the bubble flowfield. 




Comparison of the shape factor correlation of Eppler, Drela, and the present 
function. 
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The height of the bubble at transition is estimated by means of the expression 


hr = 


li 

tan 7 


(5.31) 


where tan 7 is given by Eq. (3.18). An estimate for the spreading angle is obtained 
from the experimental values reported by Birch and Eggers [1973]. In this report, 
measured spreading rates are expressed in nondimensional form as function of a 
velocity ratio defined as 


u 1 — u 2 

Ui + u 2 


(5.32) 


where U\ and u 2 are the velocities above and below a splitter plate, respectively. 
The spreading rates are normalized with respect to the maximum, which occurs at 
A = 1 . The data are all taken in shear layers developing in zero pressure gradient, 
so that their validity in the turbulent part of the bubble, where the pressure varies 
very rapidly, could be doubted. 

Making recourse to the strongly interacting nature of the flow, the above ob- 
jection can be dispelled. Specifically, this shear layer is not developing inside a 
duct with diverging wall, where the pressure gradient is imposed as a boundary 
condition of the inviscid flow. In the bubble, the amount of pressure recovered is 
strictly a function of the intensity of the turbulence: of the momentum transfer 
across the shear layer. Therefore, treating the reattaching turbulent shear layer in 
the bubble as a shear layer in zero pressure gradient with varying velocity ratios, 
with the rise in pressure a by-product with negligible feed-back, seems a reasonable 
approximation. 

In order to obtain an estimate of the spreading rate of the shear layer from the 
experimental plot given by Birch and Eggers, it should be realized that between 
transition and reattachment A varies between a value slightly greater than 1 to 1 . 
It falls, therefore, off the plot. Since the magnitude of the reverse flow is quite small 
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compared to the local edge velocity, it is not a bad approximation to assume that 
the spreading rate will be slightly larger than the maximum reported and constant. 
The following function, which depends weakly on Reynolds number, has been found 
satisfactory, 

tan 9 = .0975 + 2.5 x 10 8 T2 (5.33) 


This corresponds to a spreading angle varying from 5.71° at R = 100,000 to 8.39° 
at R = 2,000,000, measured from the parallel to the airfoil surface. Thus, the 
turbulent length of the bubble is simply obtained from 




for 

tan 9 


(5.34) 


The bubble geometry as defined by these expressions is summarized in Fig. 5-3. 

The coarseness of this approximation may seem unnecessary. In fact, it is a 
well known fact that the dividing streamline is curved downstream of transition 
and meets the airfoil surface at 90°. It is assumed here, instead, to be straight. 
In addition, the spreading of the shear layer should be measured from its bottom 
edge or from its center, the line of zero velocity, whereas here it is measured from 
the dividing streamline. This apparently wrong resolution of the bubble geometry 
is followed because the height of the dividing streamline is the only length scale 
known entirely from the upstream development. The approximation works so well 
because the extent of the turbulent part of the bubble is very short and because 
this height is an accurate characteristic length. 

Having determined ^ 2 ? Eq. (5.1) is completely specified when values are as- 
signed to the A{. The value for A\ may be linked to the length of the transi- 
tion region. In accordance with the i^-distributions measured by Horton [1967] 
as well as, more recently, by Fitzgerald and Mueller [1990], A\ is such that H 32 
grows steeply downstream of transition to a local maximum before dropping to the 
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reattachment value. This local maximum corresponds to the sharp “knee” in the 
pressure distribution which, as maintained by Russell [1978], occurs downstream 
of transition and corresponds to the location where the turbulent shear layer first 
touches the surface. The following function has been found necessary to give good 
results for leading-edge as well as mid-chord bubbles, 

Ai = 0.5 + e - 30 °( h r/c ) (5.35) 

The value of A^ is obtained by iterating on the intersection angle between the 
undershoot and the inviscid distribution. As shown in Fig. 5-4, the calculations 
between the reattachment point and this intersection are repeated with different 
values of A 2 until the angle satisfies the desired tolerance. As the drag prediction 
has been found to be quite insensitive to this tolerance, it has been relaxed to only 
1 radian to maximize computational efficiency. 

The closure relations for Co and cy, Eqs. (5.23) and (5.30), were originally 
used unchanged. There were cases, however, when the reattaching experimental 
pressure distribution could not be reproduced. This is explained as follows. Drela’s 
turbulent correlations were originally developed for separated turbulent boundary 
layers downstream of turbulent separation. The skin-friction coefficient obtained 
from these correlations in the turbulent part of the bubble is quite small in magni- 
tude, equal to or less than that in the laminar part, while the dissipation coefficient 
is very similar to that of an attached turbulent boundary layer. While such val- 
ues make sense in the slowly recirculating, constant pressure flow downstream of 
turbulent separation, they are clearly too small to reflect the effect of a turbulent 
shear layer impinging on the wall, a process for which the boundary-layer approx- 
imations are likely to break down. This claim is also supported by Navier-Stokes 
simulations, such as the one reported by Briley and McDonald [1983], which show 
a peak in negative skin friction in the turbulent part of the bubble several times its 
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value before transition, a s well as by the high peak in heat-transfer coefficient mea- 
sured at the reattachment point of bubbles developing on turbine blades [Pucher 
and Gohl, 1987]. Roberts [1980] reports a measured mean value for the dissipation 
coefficient of 0.035, which is twice the value originally proposed by Horton and also 
given by Drela’s correlations. Finally, it is not known with certainty whether or not 
the reattachment process is unsteady. Unsteadiness certainly seems likely at lower 
Reynolds numbers, as the critical value is approached, or for very large values of 
hr l c , The most effective way to capture such unsteadiness in a steady, integral 
method is to lump its effects into a higher value for Cp. 

Two additional empirical functions have been introduced. The distribution of 
skin-friction is specified by fitting a parabola to the transition point, the reattach- 
ment point (where cj — 0 by definition), and to a preassigned value, Cf min , half-way 
in between, 

c/ mi „ =-^0.0002^ (5.36) 

The higher level of dissipation is obtained by means of a multiplicative function 
to Drela’s values. This function rises quadratically from 1 at transition to a peak, 
^ Drnax ? the mean reattachment point and then exponentially decays back to 1 a 

small distance downstream of reattachment, 


'•{ 

1 + (C Drnai - 1) (^) 2 , 0 < < 1 

i + (c , D m .,-i)e 4 ^-0, *=j*>i 

(5.37) 

where the decay rate 

is 



hr 

v = 15 - 1000 — 
c 

(5.38) 

such that 

Cd = / X [ColDrela 

(5.38) 


The decay rate downstream of reattachment increases with decreasing bubble thick- 
ness at transition. This reflects the assumed slower rate of decay of the turbulence 
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(or unsteadiness) downstream of thicker bubbles. 

It is not known whether or not Roberts’s value of Cd = 0.035 is representative 
of most bubbles. Thus, although the value of Co rnax could be adjusted to match 
this experimental value, its general dependence on varying flow conditions cannot 
be inferred from it. Given the great impact of this variable on the shear layer 
development, however, it is still indispensable to obtain an estimate for at least 
the order of magnitude of its variation. Rather than offering a rigorous derivation, 
the following argument discusses flow variables that may be used to develop an 
empirical correlation between Co rnax and hj jc. 

As shown in Fig. 5-5, it is observed that, for an airfoil at a fixed a, as the 
Reynolds number decreases the bubble increases in length with the (nondimen- 
sional) edge velocity at transition remaining practically constant. The laminar sep- 
aration point calculated with the boundary-layer equations (without interaction) is 
independent of Reynolds number [Schlichting, 1979]. Since most airfoil velocity dis- 
tributions are nearly linear in the main recovery region, the amount of velocity that 
needs to be recovered in order to reach the inviscid distribution from a constant 
value at transition increases linearly as this value moves downstream. Referring 
now to Fig. 5-6, it can be seen how the main contribution to the decrease in edge 
velocity comes from the dissipation coefficient term in Eq. (5.6). Specifically, the 
decrease in velocity between transition and reattachment is mostly dependent on 
the area under the Co-distribution. As an aside, it can be seen from this figure 
that a correct modelling of c/ is not crucial. If a linear decrease in velocity with 
increasing transition length is desired, therefore, an approximately linear increase 
in this area is necessary. As the Reynolds number decreases, the ratio 

'/(*.)' t 5 - 39 ) 

which determines the separation angle, varies together with the transition length 
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Effect of decreasing Reynolds number on bubble velocity distribution 
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in such a way that the height of the bubble at transition varies almost linearly 
with t \ . This is shown in Fig. 5-7 for the bubble and the Reynolds numbers 
shown in Fig. 5-5. Finally, as the spreading angle of the turbulent shear layer is 
also nearly constant, the turbulent length of the bubble also increases linearly with 
transition length. Therefore, in order for the area under the Co-distribution to 
increase linearly over this length, Co moi should remain constant. 

Given the very approximate nature of the argument given above, it is not 
surprising that some variation in Co mox was found necessary. The function that 
has given the best agreement is 

C Dmax = 1 + ^200^ (5.40) 

This modification has led to a much greater control on the amount of pressure 
recovery between transition and reattachment such that any experimental pressure 
distribution can now be reproduced simultaneously with the correct growth in 62 . 
This function is shown in Fig. 5-8 together with the variations of A\ and c/ mtn . 

Model Flowchart 

Having described each part of the bubble separately, the scheme used by the 
present method of predicting the development of laminar separation bubbles is now 
summarized by the flow diagram shown in Fig. 5 - 9 . Starting with the inviscid veloc- 
ity distribution over an airfoil, the bubble model is invoked when laminar separation 
is predicted. After removal of the Goldstein singularity, ( Rs 2 )s is determined and, 
based on the inviscid velocity gradient at the laminar separation point, an initial 
estimate of Gaster’s parameter, P, is made. These two parameters are necessary 
to estimate the angle that the separating streamline makes with the surface. The 
velocity distribution in the plateau region is prescribed by means of the velocity 
plateau function which depends both on P and on the matching of its slope to that 
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of the inviscid velocity distribution at the laminar separation point. The separated 
shear layer development can thus be calculated in the direct mode using the mo- 
mentum and energy integral equations. The laminar length of the bubble extends 
to the point where the amplification factor n = 9. This length times the tangent 
of the separation angle gives the thickness of the bubble at transition. From this 
height, knowledge of the spreading angle of the turbulent shear layer allows the 
calculation of the turbulent length of the bubble. The value of the inviscid ve- 
locity corresponding to this known reattachment location can be used to obtain a 
new value for P , such that the laminar calculations can be iterated until P reaches 
a fixed value. Upon convergence, the shear layer development in the turbulent 
part of the bubble is calculated by prescribing the distribution of iJ 32 and solving 
the integral boundary-layer equations in the inverse mode together with turbulent 
closure relations based on Drela’s but modified to model better the reattachment 
process. Upon the intersection of the undershoot with the inviscid distribution, 
the boundary-layer development is calculated in the direct mode using the inviscid 
pressure distribution to drive Drela’s unmodified non-equilibrium turbulent to the 
trailing edge, where the drag is obtained with the Squire- Young [1937] formula. 

The equations used in the different parts of the bubble model are summarized 
in the Appendix. 
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Chapter 6 
RESULTS 


In this chapter, the bubble model will be tested by comparing its predictions to 
available experimental measurements. These include mostly drag polars and pres- 
sure distributions for several airfoils but also two sets of Laser-Doppler Velocimetry 
measurements inside the bubble. The Reynolds number of these tests ranges from 
2,000,000 down to 100,000. 


NACA 66,-018 Airfoil 

The NACA G63-OIS airfoil was one of the first to be tested for which the effect 
of the bubble on the pressure distribution could be clearly seen. Figs. 6-1 and 6-2 
show a comparison between the predicted pressure distribution and boundary-layer 
developments and those measured by Gault [ 1955 ]. It can be seen that away from 
the bubble the inviscid pressure distribution is quite satisfactory. To the right of the 
plot is a blow-up of the bubble pressure distribution. The two asterisks represent 
the viscous separation and the reattachment points. The local inverse solution 
near separation with H12 prescribed employed to remove the separation singularity 
results in a rounding of the discontinuity in the velocity gradient upstream of sep- 
aration and in some upstream influence of the bubble on the pressure distribution. 
Although this is achieved by purely numerical means, the correct local behavior 
seems well captured. The growth of the amplification factor is plotted along the 
airfoil surface itself in units of percentage chord. This curve is plotted as function 
of xjc rather than arc-length for consistency with the pressure distribution. Thus, 
the end of the curve, at n = 9 , lies a distance equal to 9 % of the airfoil chord above 
the airfoil s y-coordinate corresponding to the transition location. For this case, all 
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of the amplification occurs after laminar separation. 

Fig. G-l(b) contains all the relevant inputs and outputs of the boundary-layer 
analysis. The inviscid velocity distribution is provided for reference together with 
the Reynolds number and angle of attack. The calculated bubble is drawn on 
this plot. The summary of the viscous analysis as printed out by the Eppler and 
Somers program is printed to the right, below the label indicating that a bubble 
analysis has been performed with transition at n = 9. This includes the extents 
of turbulent and separated arc lengths normalized with respect to the chord, as 
well as the drag coefficients from the upper and lower surfaces. The lift and the 
total drag coefficients are printed next, and the bubble lengths on the upper and 
lower surfaces come last. The Eppler boundary-layer development plot is shown 
below. This plot is especially useful during the design process. Finally, all five 
boundary-layer variables are plotted as functions of x/c, with the asterisk denoting 
the reattachment point. The removal of the Goldstein singularity can be clearly 
seen in the smooth developments through the separation point. 

This figure serves to illustrate several points. Since the airfoil is symmetrical 
and a — 0, by suppressing the bubble model on the lower surface the present 
prediction using the model and Drela’s turbulent boundary-layer method can be 
compared with the original Eppler turbulent boundary-layer analysis starting at 
the laminar separation point. With the exception of the bubble region, the two 
methods give very similar results, as expected. The difference in separation point 
locations with and without Goldstein’s singularity can be seen from the plot of 
H 12 vs. x/c. At the bottom left of the figure Eppler’s boundary-layer development 
plot shows that laminar separation for this airfoil occurs at the boundary with 
natural transition according to the modified transition criterion Eq. (4.5). The e n 
analysis, however, indicates that n at separation is still quite small. This apparent 
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inconsistency of Eppler’s transition criterion may be resolved by realizing that it was 
calibrated using mainly drag data. For this case, the drag of the upper and lower 
surfaces is identical. Therefore, the bubble does not seem to cause any deterioration 
in airfoil performance. The experimental values of momentum thickness inside the 
bubble, derived from Gault’s data by Roberts [1980], support this interpretation. 
As shown more clearly in Fig. 6-2, a bubble may not necessarily lead to an increase 
in momentum thickness greater than if transition had been assumed at laminar 
separation. In fact, as will be shown later, below a certain length a mid-chord 
bubble appears to reduce the airfoil drag. The growth in momentum thickness in 
the laminar part of the bubble is clearly evident. The calculated transition point is 
a few percent chord too far downstream. As the length of the plateau in the laminar 
part shown in Fig. 6- 1(a) is quite close to the experimental, transition may indeed 
start before any change in the pressure is observed, as maintained by Russell [1978]. 

NASA NLF(1)-1015 Airfoil 

This airfoil was recently tested in the NASA Langley Low- Turbulence Pressure 
Tunnel. Drag polars calculated from force measurements for lift, wake surveys 
for drag, and also detailed pressure distribution measurements were obtained at 
R= 2, 000, 000, 1, 000, 000, 700, 000, and 500, 000. The profile was designed for use 
on a high-altitude long-endurance RPV and is therefore characterized by a large 
aft-loading to achieve the high ct requirement. It was also designed to minimize the 
effects of the bubbles at the design conditions of R = 2,000,000 for a high-speed 
dash (bottom of the low-drag bucket) and R = 700,000 for maximum endurance 
(top of low-drag bucket). The present predictions are compared with the R = 
500,000 data since here the effects of the bubble are most clearly seen. 

Fig. 6-3 shows the aerodynamic characteristics of this airfoil on the plot as gen- 
erated by the original Eppler and Somers program. Shown is the original prediction 
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Fig. 6-3 Aerodynamic characteristics of the NASA NLF(l)-1015 airfoil obtained with the 
original Eppler and Somers program compared with those obtained with the present 
bubble model and with the interactive program XFOIL of Drela, R — 500, 000. Data 
are from the Low-Turbulence Pressure Tunnel, NASA LaRC, 1987. 
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obtained by assuming transition at laminar separation, the prediction of the present 
model, Drela’s XFOIL results, and the experimental data from LTPT. Part of the 
difference in drag prediction between the original program and the present version 
is due to the different turbulent boundary-layer methods employed, as shown in Fig. 
2-4. In fact, comparing Figs. 2-4 and 6-3, it can be inferred how at the top of the 
bucket the bubble leads to a drag reduction over what is calculated by assuming 
transition at laminar separation. In any case, the present formulation is able to re- 
produce the measured data with excellent accuracy. At the upper and lower limits 
of the polar, the onset of strong global interaction cannot be neglected. At these 
conditions, the present lift and drag predictions are poorer. The plot also contains 
the lift and moment curves as functions of a as well as the transition and turbulent 
separation locations as functions of x/c. In actuality, although the program labels 
the axis as “x/c,” the independent variable is really “1 — St ur b/c ,” a rather more 
cumbersome variable. Given that there is usually very little difference between the 
two variables, x/c is used for ease of presentation. Since the original version of 
the program assumes transition at the laminar separation point, the difference be- 
tween the transition curves for the two analyses represents the laminar length of the 
bubble. It can be seen how the bubble decreases in size and eventually disappears 
as the pressure distribution upstream of laminar separation becomes increasingly 
adverse. This will be clearly shown in subsequent plots. Because the XFOIL tran- 
sition locations are given in terms of actual x/c, they appear to occur downstream 
of transition as calculated by the present model. As will be shown below, however, 
they occur upstream. 

Matching the experimental drag polar does not by itself guarantee an accurate 
prediction. The pressure distribution and the boundary-layer development should 
also be compared to experimental data. Unfortunately, it is much more difficult 
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to measure these quantities. Consequently, few data sets are available. For the 
NLF(1)-1015 airfoil, the details of the bubble flowfield can be checked only through 
the pressure distribution, since no boundary-layer data were taken. Figs. 6-4 to 
6-6 show typical comparisons with the measured pressure distributions together 
with the corresponding calculated developments. Fig. 6-4(a) shows the pressure 
distribution corresponding to the lowest value of ci on the polar together with the 
^^OIL prediction. Two very different bubbles can be seen, one at the mid-chord 
on the upper surface and the other at the leading edge on the lower surface. Both 
bubbles are well approximated by the model. XFOIL gives a slightly shorter bubble 
on the upper surface and a very slight perturbation on the pressure distribution on 
the lower. Fig. 6-4(b) shows the boundary-layer development. It is interesting 
to see how large the values of the shape factors can become inside leading-edge 
bubbles. These, in any case, do not seem to contribute much to increasing the drag 
of the airfoil. The shape of the sin(l/x) function can be clearly recognized in the 
mid-chord bubble. The values for C& rnaT are very similar, the greater increase in 62 
of the upper-surface bubble being largely due to the longer extent of its turbulent 
part. This bubble, in fact, is much thicker at the transition point than the leading- 
edge bubble. The distribution of Cf seems plausible. In any case, in this region this 
variable has little impact on the boundary-layer development or on the pressure 
distribution. 

Fig. 6-5(a) shows the pressure distribution at an angle of attack corresponding 
to the middle of the airfoil polar. In an attempt at matching the pressure gradient 
along the bubble, the inviscid angle of attack is one degree less than the experi- 
mental. It can be seen, however, that the strong trailing-edge interaction induced 
by the large aft-loading prevents the matching of the gradients on both surfaces 
simultaneously. In any case, the model reproduces the measured pressure distribu- 
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tion remarkably well. Both the upper and lower surface bubbles are quite long and 
may be expected to be thick. In fact, as can be seen in Fig. 6-3, at this condi- 
tion the drag due to the bubble is highest. The bubbles predicted by XFOIL are 
a little short although the same value of n was used. At this angle of attack, n is 
starting to be amplified upstream of laminar separation on the upper surface. The 
boundary-layer development is shown in Fig. 6-5(b), where the step in 82 in the 
turbulent part of the bubble is evident. 

Fig. 6-6(a), finally, corresponds to the top of the bucket, where the upper 
surface bubble is about to disappear and the lower surface contributes little to 
the total drag. The short transition length in this case is believed to be a direct 
consequence of the destabilizing effect of the adverse pressure distribution upstream 
of separation. This is clearly shown by the distribution of n along the upper surface. 
In Fig. 6-6(b), the viscous analysis summary shows how the lower-surface drag is 
only a fifth of the upper-surface drag in spite of a 16%c long bubble. It is particularly 
interesting to observe how the upper-surface bubble is shrinking while preserving 
its proportions. This indicates that the correct scaling for the bubble has been 
identified. At a slightly higher angle of attack, the transition point corresponds to 
the laminar separation point. Beyond this condition, the transition point precedes 
the separation point and travels upstream until the rise of the suction peak again 
leads to laminar separation before transition and to the formation of a leading-edge 
bubble. As far as the effects on transition are concerned, therefore, it appears that 
the a destabilizing pressure distribution is entirely analogous to a rise in Reynolds 
number. 


Eppler E387 Airfoil 

This airfoil was designed more than twenty years ago and is intended for use 
on model gliders. It also was recently tested in the NASA Langley Low- Turbulence 
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model and with the interactive program XFOIL of Drela, R = 300,000. Data are 
from McGhee et al. [1988]. 
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Pressure Tunnel [McGhee et al., 1988]. Figs. 6-7 to 6-10 show comparisons between 
measurements and calculations at R — 300,000. Fig. 6-7 shows the aerodynamic 
characteristics. While the original program underpredicts the drag by 10 or 20 
counts, the new version is quite accurate. Also shown is the prediction with Drela’s 
XFOIL program. Even though XFOIL was run with n = 12, the drag is still 
underpredicted. As the bubble in this case is as long, if not longer, than the exper- 
imental, the small drag values can only be a consequence of a too small value for 
the dissipation coefficient in the turbulent part of the bubble. 

Fig. 6-8(a) shows the pressure distribution for the lowest point on the polar. 
The slight drag overprediction at low ct in Fig. 6-7 is caused by the too steep 
bubble recovery that can be seen on the upper surface. Fig. 6-8(b) shows the 
boundary-layer development. Fig. 6-9 shows the results at a = 1.5°. Fig. 6- 
10(a) shows a leading-edge bubble at a high c*. It is quite similar in shape to the 
slight perturbation in the experimental pressure distribution, including the small 
undershoot. Better agreement might be obtained with slightly different profiles, for 
instance the Green profiles. In fact, the early transition is a consequence of the very 
high values for the shape factor, shown in Fig. 6-10(b), which seem unrealistic. 

Fig. 6-11 shows the drag polar at R = 200,000. Again, the prediction is 
excellent. Figs. 6-12-6-14 are some characteristic analyses. Fig. 6-15 is a limiting 
case, very near the critical Reynolds number. Below R = 100,000, the model 
breaks down due to a transition length which extends beyond the trailing edge. 
At this Reynolds number the prediction is not as good, as shown in Figs. 6-16 
to 6-20. The pressure distribution shown in Fig. 6-16(a) highlights some of the 
limitations of the model. The error due to the envelope method increases as the 
Reynolds number decreases. Thus it could explain the disagreement between the 
measured and predicted transition predictions for the upper surface bubble. The 
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Fig. 6-12 (a) Comparison of predicted and measured pressure distribution for the Eppler E387 
airfoil. Experimental a = —2°. Data are from McGhee et al. [1988], 
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Fig. 6-14 (a) Comparison of predicted and measured pressure distribution for the Eppler E387 
airfoil. Experimental a = 7°. Data are from McGhee et al/ [1988]. 
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Fig. 6-15 Aerodynamic characteristics of the Eppler E387 airfoil obtained with the original 
Eppler and Somers program compared with those obtained with the present bubble 
model and with the interactive program XFOIL of Drela, R = 100,000. Data are 
from McGhee et al. [1988]. 
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Fig. 6-16 (a) Comparison of predicted and measured pressure distribution for the Eppler E387 
airfoil. Experimental a = -2°. Data are from McGhee et al. [1988], 
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Fig. 6-17 (a) Comparison of predicted and measured pressure distribution for the Eppler E387 
airfoil. Experimental a = 0°. Data are from McGhee et al. [1988]. 
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Fig. 6-18 (a) Comparison of predicted and measured pressure distribution for the Eppler E387 
airfoil, a = —5°. Dot-dashed line is XFOIL prediction at a = 2°. Data are from 
McGhee et al. [1988]. 
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Fig. 6-19 Cont. (b) Viscous analysis summary and boundary-layer development for the Eppler E387 
airfoil, a = 4°. 
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Fig. 6-20 (a) Comparison of predicted and measured pressure distribution for the Eppler E387 
airfoil. Experimental a = 7°. Data are from McGhee et al. [1988]. 
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lower surface bubble is of a very difficult type for the model. In fact, the velocity 
distribution along the bubble deviates significantly from a straight line, such that 
the height at the transition point is probably overestimated by Eq. (5.31). As a 
consequence, the predicted turbulent length is also too large leading to a too great 
pressure recovery. A higher-order dependence of h j on the curvature of the velocity 
distribution should be developed. As discussed by Gaster, such deviation from a 
linear recovery may be well represented by a P 2 -term in the correlations. As the 
pressure distribution upstream of laminar separation becomes increasingly adverse, 
as shown in Figs. 6-17(a)-6-20(a), the predicted transition point moves downstream 
relative to the experimental. This is consistent with the conclusions of Chapter 4 
about Drela’s e n method. In Fig. 6- 18(a), the effect of the interaction on the 
growth of n can be clearly seen in the XFOIL analysis. At this Reynolds number, 
in fact, the bubble is so large that its effect on the pressure distribution upstream 
of laminar separation cannot be neglected anymore. The steeper growth of n that 
results is responsible for the early transition. This effect is probably independent 
of the e n method employed such that a correlation between n and R should be 
developed for use with interactive methods. In Figs. 6-16(b)-6-20(b) the details 
of the boundary-layer development inside the bubble are most clearly evident. To 
conclude the comparisons, the last two airfoils to be discussed, Figs. 6-21 and 6- 
22, were tested for pressure distributions and boundary-layer developments using 
Laser-Doppler Velocimeters. 


NACA 65-213 Airfoil 

An NACA 65-213 airfoil was tested by Hoheisel et al. [1984] at DLR, Ger- 
many. Pressure coefficients and boundary-layer developments obtained with a 
Laser-Doppler Velocimeter are given for the upper surface at zero angle of attack 
and R = 240,000. Fig. 6-21(a) shows the airfoil and the pressure distribution. Al- 
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Fig. 6-22 (a) Comparison of predicted and measured pressure distribution for the Wortmann FX 
63-137 airfoil. Experimental a = 7°. Data are from Brendel and Mueller [1986], 
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Fig. 6-22 Cont. (b) Viscous analysis summary and upper surface boundary-layer development for the 
Wortmann FX 63-137 airfoil, a = 5° compared to the Laser Doppler Velocimeter 
measurements of Brendel and Mueller [1986], a = 7°. 
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though the inviscid analysis is performed at one degree less than the experimental, 
the pressures could not be matched since the coordinates used were not the exper- 
imental ones and, moreover, were generated by scaling those of an NACA 65-210 
airfoil. As explained in Abbot and Von Doenhoff [1959], this is an approximate 
procedure which can be used in place of the exact method of Theodorsen if the 
change in thickness is small. Finally, the turbulence intensity of the wind tunnel 
used is TiZoo = 0.2%. This is to be contrasted with that of LTPT, which is less than 
0.02%. Both the shape of the bubble distribution and the transition location appear 
to be affected by such a high freestream turbulence intensity. In fact, the bubble 
recovers a significant amount of pressure in the laminar part, followed by a fairly 
gradual steepening into the fully turbulent recovery. As shown in Fig. 6-21(b), the 
experimental transition point occurs about 5%c upstream of the predicted. In any 
case, the shape factor and momentum thickness developments are reproduced very 
well. 


Wortmann FX 63-137 Airfoil 

The last of the comparisons is a prediction of the pressure distribution and 
bound ary- layer development over a Wortmann FX 63-137 airfoil at Ft = 100, 000. 
The data were taken by Brendel and Mueller [1988] with an LDV at a - 7 . 
Fig. 6-22(a) shows the pressure distribution and Fig. 6-22(b) the upper surface 
boundary-layer development. Although the calculated trailing-edge value of mo- 
mentum thickness is 0.028, the plot is cut off as shown in order to view more clearly 
the bubble region. The comparison is very good. 
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Chapter 7 
CONCLUSION 


Summary and Conclusions 

The correct scaling parameters for the bubble have been determined: P, 

(R$ 2 )s, and hr / c. The generality of the model relies on having understood and 
correctly approximated the dominant physical processes in the bubble flowfield. 
Having identified what the bubble depends on, parameters characteristic of this 
particular problem have been utilized to construct an algorithm that conforms to 
the flow development itself. Thus, the strong interaction leading to the pressure 
plateau in the laminar part is well represented by DU(P ), a relationship that re- 
flects the dependent role of the pressure recovery. The e" method coupled with 
Wortmann’s correlation for the separation angle as a function of P and (R$ 2 )s al- 
lows the transition location and the height of the bubble at transition to be found. 
^ T / c t in. turn, determines the length of the turbulent part such that even a rough 
estimate of the turbulence intensity in the reattachment region allows an accurate 
prediction of the pressure distribution and the momentum thickness growth. 

The feasibility of a semi-empirical approach has finally been established. After 
fifty years of unsuccessful attempts at developing a general semi-empirical bubble 
model, the identification of the correct scaling parameters for the bubble has finally 
made such an approach possible. This has brought two main advantages: a method 
more efficient than any interactive boundary-layer of finite-difference algorithm and 
a deeper understanding of the bubble flowfield without which, in fact, the method 
could not have been developed. 

A general, accurate, and computationally efficient laminar separation bubble 
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model has been developed. Although very simple, the model is able to reproduce 
the effects of vastly different bubbles over the whole Reynolds number range in 
which bubbles form. On a VAXstation 3100, the original version of the Eppler and 
Somers program takes approximately ten seconds to generate a drag polar, which 
typically is defined by fourteen angles of attack. The program with the bubble 
model takes approximately one minute for the same analysis. Drela’s XFOIL takes 
approximately thirty minutes, partly because smaller steps in angle of attack must 
be taken to ensure convergence. 

Suggestions for Future Work 

The simplicity of the model brings a few drawbacks. While its response to 
large changes in the controlling flow conditions is correct to lowest order, its sen- 
sitivity to slight variations in pressure distribution is limited. More importantly, 
its performance near limiting conditions such as high angles of attack, very low 
Reynolds numbers, or very unusual pressure distributions is not reliable. These 
deficiencies can be traced to the “stiffness” of the model: by effectively integrating 
parts of the flowfield in an approximate way, the pointwise flexibility of the original 
governing equations is lost. Whereas it appears that the bubble flowfield itself has 
been properly “integrated,” its effect on the airfoil pressure distribution has not. 
An interactive method should therefore be incorporated in the Eppler and Somers 
program in order to be able to analyze those limiting cases where the present model 
is likely to fail. 

The principles of conservation of mechanical energy and momentum that are 
invoked to justify the pressure plateau in the laminar part might be used to deduce 
what the value of G should be in the Green profiles. Since the value of skin friction 
appears to have such a small impact on the bubble development, the correlations 
based on these profiles may lead to a more accurate calculation of the separated 
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laminar shear layer. 

Another concern for future work on the improvement of the model is the 
transition prediction method. The method of Stock and Degenhart should be 
implemented in place of Drela’s but its accuracy should nonetheless be tested 
against “exact” stability analyses of boundary-layer developments calculated by 
finite-difference methods. The influence of the interaction on the growth of n de- 
serves particular attention. It appears from the XFOIL analyses that the correct 
transition location can be matched by a value of n that increases with decreasing 
Reynolds number. It should therefore be possible to develop an empirical function 
n(R) that reflects this trend. 

Finally, very detailed measurements of the reattachment process are necessary 
to understand it better and to provide a better estimate for Co in the turbulent 
part of the bubble. Specifically, the mean flow and the turbulent stresses should be 
measured in order to deduce the actual distribution of dissipation coefficient in the 
turbulent part of the bubble. 


APPENDIX 


MODEL SUMMARY 


In this appendix, all the equations used in the present version of the bubble 
model are summarized. This will facilitate both the understanding and the repro- 
duction of the model. 


Laminar Boundary-Layer 

The inviscid velocity distribution is used to drive the laminar boundary-layer 
development, which is found by integrating the following system of equations. 


Governing Equations: 


d 6 2 
ds 


c f (ff + 2 /*^ 
j-(H 12+ 2 — 


dS 3 _ J 3 dU 

ds C ° 3 U ds 


(A.l) 

(A.2) 


Closure correlations: 


77 12 


(25.7157857477 32 - 89.58214201)77 32 
+79.87084472, 77 3 2 > 1.7258 

\Aff32 - 1.515095[(— 227. 18220 77 32 
+724.55916)77 32 - 583.60182], H 32 < 1.7258 


(A3) 



(2.22168722297732 - 4.226252829)77 32 

+1.3723907030, 77 32 > 1.7258 

[(— 0.0317285065577i2 + 0.3915405523)77 12 

-1. 686094798] 77 12 + 2.512588652, 77 32 < 1.7258 


Rs 2 C d = (6.837796177 32 - 20.521 103)77 32 + 15.707952 


(A5) 


Laminar Separation: 


77 32 = 1.515095 


(A6) 
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Laminar Part of the Bubble 

The development of the shear layer in the laminar part of the bubble is cal- 
culated by integrating Eqs. (A.l) and (A. 2) along with the following expressions. 
The removal of the Goldstein singularity at laminar separation is discussed at the 
end of Chapter 3. 


Separation angle: 


where 


tan 7 — 


64 P 


( R 6*)i 


(Rs 2 )s = R 


Us {62 )s 
U 00 c 


P = R 




1 2 


MU/Uoo) 

A(s/c) 


Laminar pressure recovery 
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DU U s 
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]} 


where 


DU 


= f 0.0610 

1 0.0152 
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Closure correlations: 
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Transition 


Amplification factor: 


n(^) = [ 

J SQ 


dn 

dR $ , 


(flu) 


m(i?i 2 ) + 1 l(#12) 
2 <5 2 (-s) 




where 


dn 

dRs„ 


(H 12) 


P») 

m(/f 12 ) 


= 0.01[{2.4F 12 - 3.7 

+ 2.5 tanh[l. 5(^,2 - 3.1)]} 2 +0.25]* 

6.54tf 12 - 14.07 
I 72 

n 12 


0.058(fTi 2 - 4) 2 1 
(H n - 1) - 0.068 i 


•so is defined as the location where R$ 2 = Rs 2o , where 


logi 0 [i? i20 (tfi2)] = 


+ 


1.415 
H 12-1 
3.295 


- 0.489 


tanh 


20 


H 12 -1 


- 12.9 


# 12-1 


+ 0.440 


Laminar length of the bubble: 


£i@n - 9 


Bubble height at transition: 


hr — 


h 

tan 7 


Turbulent Part of the Bubble 
Spreading angle of turbulent shear layer: 


(A. 15) 

(-4-16) 

(-4.17) 

(-4.18) 

(-4.19) 

(A. 20) 
(-4.21) 


tan# = .0975 + 2.5 x Kr’R 


(A. 22) 
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Turbulent length of the bubble: 


?2 = 


tan 6 


Shape factor distribution: 


where 


#32 (y) = sin 


7T 

L V 


#32 ~ (#32)tt 


+A, [(# 32 )r - (#32)*] 


#32 ~ 


y = l 3 - Vo ) ° + Vo 


Vo - 


7T 


37t — sin 1 ( 1 /^4 1 — 1) 

f ( 5 ~ st) ft 2 s < s-r 

1 [(^ — sr )/^2 — 1]5F + 1 8>Sv. 

SF — •y/ Ai /A 2 

A x = 0.5 + e- 300 ^ r / c ) 


(A. 23) 


(A. 24) 


(A. 25) 

(A. 26) 
(A. 27) 

(A. 28) 
(A.29) 
(A. 30) 


and A 2 is iterated upon until the undershoot merges with the inviscid velocity 
distribution. 

The above shape factor distribution is used to drive the integral boundary-layer 
equations in the inverse mode. 


Governing equations: 

dU 

ds 

d 6 2 
ds 


(c f /2)H 32 -C D +6 2 


d# 32 

ds 


U 


^ 2 # 3 2 (#12 — 1 ) 
#12 + 2 
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(A.31) 
(A. 32) 
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C T ds 


(A. 33) 
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where 
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Closure correlations: 
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e/_,„ = -^0.0002^ (A.48) 

The c/ distribution is obtained by fitting a parabola through c/ T , c/ mfn , and c/ = 0 
at ^2 ■ 


Turbulent Boundary Laver 

The development of the turbulent boundary layer downstream of the intersec- 
tion of the undershoot with the inviscid velocity distribution is calculated in the 
direct mode by prescribing the inviscid velocity distribution to drive Eqs. (A.l), 
(A. 2), and (A. 33). This system of equations is complemented by closure correlations 
(A. 34), (A. 35), and (A.40)-(A.47). The skin-friction coefficient is obtained from 


0.3e” 1,33 ^ 12 

c f = + 0.00011 

(log^,) 1 - 74+on " 1 ’ 

?urbulent separation: 

H 3 2 = 1.505 + 


tanh ( 

4- 

V 

0.875 ) 


400 


(A. 49) 


(A. 50) 


Squire- Young formula: 


Drag Calculation 





2 . 5 + 0.5 H 13te 


Cd = 


C 


(A.51) 
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